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ABSTRACT 

We present measurements of the auto- and cross-frequency power spectra of the cosmic infrared 
background (CIB) at 250, 350, and 500 |J.m (1200, 860, and 600 GHz) from observations totaling 
~ 70 deg 2 made with the SPIRE instrument aboard the Herschel Space Observatory. We measure a 
fractional anisotropy SI /I — 14 ± 4%, detecting signatures arising from the clustering of dusty star- 
forming galaxies in both the linear (2-halo) and non-linear (1-halo) regimes to unprecedented levels; 
and that the transition from the 2- to 1-halo terms, below which power originates predominantly 
from multiple galaxies within dark matter halos, occurs at kg ~ O.larcmin -1 (I ~ 2200). New to 
this paper is clear evidence of a dependence of the Poisson and 1-halo power on the flux-cut level of 
masked sources — suggesting that some fraction of the more luminous sources occupy more massive 
halos as satellites. We measure the cross-correlation power spectra between bands, finding that bands 
which are farthest apart are the least correlated, as well as hints of a reduction in the correlation 
between bands when resolved sources are more aggressively masked. In the second part of the paper 
we attempt to interpret the measurements using a novel alternative to the standard halo model. With 
the aim of fitting simultaneously with one model the power spectra and number counts in all bands, 
we find that this is achievable by invoking a luminosity-mass relationship, such that the luminosity- 
to-mass ratio peaks at a particular halo mass scale and declines towards lower and higher mass halos. 
Our best-fit model finds that the halo mass which is most efficient at hosting star formation in the 
redshift range of peak star-forming activity, z ~ 1 — 3, is log(Af pea k/Mo) ~ 13.9 ± 0.6, and that the 
minimum halo mass to host infrared galaxies is log(M m i n /M Q ) ~ 10.8 ± 0.6. Despite the clear step 
forward made with this model, we find that it does not fit the data to full statistical significance, and 
we suggest physically motivated modifications for future models. 

Subject headings: cosmology: cosmic microwave background, cosmology: cosmology: observations, 
submillimeter: galaxies - infrared: galaxies - galaxies: evolution - cosmology: 
large-scale structure of universe 
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1. INTRODUCTION 

Star formation is well traced by dust, which ab- 
sorbs the UV/optical light produced by young stars 
in actively star-forming regions and re-emits the en- 
ergy in the far-infrared / submillimeter (FIR/submm; e.g., 
ISavage fc Mathislll979l ) . Roughly half of all starlight ever 
produced has been re processed by dusty star-forming 
galax ies (DSFGs; e.g.. lHauser fc Dwck|[200lt iDole et at] 
l2006l l. and this emission is responsible f or the ubiq- 
uitou s cosmic infrared b ackground (CIB; iPuget et alj 
119961 iFixscn et all I1998D . The mechanisms responsi- 
ble for the presence or absence of star formation are 
highly dependent on the s pecifics of the local env iron- 
ment (e.g., major mergers: iNaravanan et al.l 120101 con- 
densation or cold accretion: iDckel et al.ll2009l photoion- 
ization heating , supernovae, actiye gala c tic nuclei, and 
virial shocks: iBirnboim &: Dckcl 120031 : iGranato et alj 
120041 iBowcr ct all I2006D . ancHis a result, DSFGs are 
likely to be non-trival tracers of the underlying dark mat- 
ter distribution. In turn, the specifics of the galaxy dis- 
tribution, which can be determined statistically to high 
precision by measuring their clustering properties, are 
valuable for models of galaxy formation. However, mea- 
suring the clustering of DSFGs has historically proven 
difficult to do. 

Owing to the relatively large point spread functions 
(PSF's) of ground-, balloon-, and space-based submil- 
limeter observatories, coupled with very steep source 
counts, maps at these wavelengths are dominated by con- 
fusion noise. For the 250 |J.m channel on Herschel, for 
example, this means that no matter how deeply you ob- 
serve a field, without some sort of spatial deconvolution 
at best only ~ 15% of the flux de nsity will be resolve d 
into individually detected galaxies ((Oliver et alj|2010bl) . 
Add to that the fact that the r cdshift distribution of 
DSFGs is relatively br oad (e.g., IChapman et al.1 [20051 : 
iBethermin et al.ll2012bl . Casey et al. in prep.), cluster- 
ing measurements of resol ved sources have consequently 
had limited success (e .g., iBlain et al.l 120041 IScott et al.1 
120061 : IWcifi et "al]|2009l) , and somewhat contradicto ry re- 
sults (e.g.. lCoorav et al.ll2010UMaddox et al.ll2010H . 

The remaining intensities in the maps appear as 
fluctuations, or anisotropies, in the CIB. Contained 
in CIB anisotropies (or CIBA) is the clustering pat- 
tern, integrated over luminosity and redshift, of all 
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DSFGs — including those too faint to be resolved. 
And analogous to the two-point function typically 
used to estimate the clustering of resolved galaxies, 
the power spectrum of these intensity fluctuations is 
a probe of the clustering properties of those galax- 
ies (e.g.. iBond fe EfstathTo^fl98l IScott fe Whitl ll99l 



iKnox et all 120011: iNegrcllo et al.M2007D. Initial power 



spect r um measurements fro m Spitzer ( Grossan fe SmootJ 
20071 lLagachc et al.l 12001. BLAST (IVicro et all 120091 



Haiian et al.1 I2012D. A CT (jDunklev et all |2011|) . and 
SPT (|Hall et all |2010D found a signal in excess of 



Poisson noise originating from the clustering of DS- 
FGs, but were limited to measuring the galaxy bias 
in the linear regime, rather than their distribution 
within dark matter h alos. Subsequent m easurements 
from Herschel /SPIRE (lAmblard et al.ll201ll) . and Planck 
(|Planck Collaboration et alj|2011bl ) were able to isolate 
the linear and non-linear clustering signals, but the two 
groups found that their measurements agreed only after 
correcting for multiple systematics. 

Power spectra can be interpreted with modeling 
frameworks in much the same way as is done for 
two-point function measurements of resolved sources. 
Among the most commo nly adopted models are so- 
called "halo models" fe.g.. IScliakll200(l ICoorav fe Shethl 
EOO I which use halo occupa t ion distributions (HODs ; 
e.g., iPeacock fc Smith! 120001 : IScoccimarro et al.l 1200 il l 
to statistically assign galaxies to dark matter ha- 
los in order to re-create observed clustering measure- 
ments. Halo models have b een adopted to in terpret 
CIBA spectra from BLAST (I Viero "et al.1 120091) Her- 



Xia et al.ll2012l). and Planck dPlanck Col 


aboration et al. 


2011bl Pcnin et al.ll2012al IShang et al.l 


20121 IXia et al. 



Precisely measuring the CIBA power spectra and 
decoding the information contained within them is a 
rapidly growing field, and it is also the focus of this pa- 
per. First and foremost, we aim to advance the field 
by providing state-of-the-art measurements of the auto- 
and cross-frequency power spectra of CIB anisotropics 
at 250, 350, and 500 urn, spanning angular scales 0.01 < 
kg < 2 arcmin" 1 (or 350 < £ < 45,000) (§ HJ. With the 

addition of more than four t imes the area, we extend 
the efforts of lAmblard et all (120111 ) — who definitively 
resolved a signature of non-linear clustering on small 
scales — by illustrating how the strength of the non- 
linear clustering signal depends strongly on the flux-cut 
level of masked s ources (§ 14.4.11). W e improve on the 
efforts of BLAST (jHaiian et al.ll2012H by measuring the 
cross-frequency power spectra and estimate the level of 
correlation between bands (§ 14. 5|) . 

We then attempt to interpret our measurements with 
a halo model whose main feature is to tie the luminosi- 
ties of sources to their host halo masses (§ 15. 1|) . Our 
model fits the power spectra in each band and measured 
P(D) number counts of sources simultaneously, thereby 
introducing a new level of sophistication to the body of 
existing halo models in the literature. When required, we 
adopt the concordance model, a flat ACDM cosmology 
with n M = 0.2 74, ft A = 0.726, H = 70.5 km s" 1 Mpc" 1 , 
and cr 8 = 0.81 (jKomatsu et al.ll201lH . 
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TABLE 1 

Map Properties of the HerMES Fields 



Field Name 


Area 


lcr Noise 


Repeats 


Scan Speed 




(deg 2 ) 


(MJy sr" 1 ) 


No. 




BOOTES 


11.3 


1.11, 0.61, 0.29 


9 


Parallel 


CDFS-SWIRE 


12.2 


1.00,0.57,0.27 


5/20 


Parallel/Fast 


ELAIS-Sl 


8.6 


1.12, 0.63, 0.31 


8 


Parallel 


LOCKMAN-SWIRE 


15.2 


1.08,0.60,0.29 


4/20 


Parallel/Fast 


XMM-LSS 


21.6 


1.15,0.68,0.52 


8 


Parallel 



Note. — Total noise {la including confusion) are given at 250, 350, 
and 500 u-m, respectively. Repeats arc defined as the number of times a 
field has been observed in two orthogonal passes, thus one repeat equals 
two passes. CDFS-SWIRE and LOCKMAN-SWIRE were observed partially in 
Parallel mode, and partly in Fast mode. The scan speed of the telescope 
is either 20arcsecs _1 (Parallel) or 60arcsecs~ 1 (Fast). 

2. DATA 

The primary data set for this work comes from the 
Herschel Multi-tie red Extragalactic Survey 1 (HerMES; 
lOliver et al.ll2012T ). a guaranteed ti me (GT) key projec t 
of the Herschel Space Observatory (|Pilbratt et a,l.|]2010f > . 
We use sub millimeter maps o bserved with the SPIRE 
instrument (jGriffin ct al.ll2010D at 250 350, and 500 \xm. 
We a lso use reprocessed 100 urn IRAS (jNeugebauer et akl 
Il984l ) maps in order to quantify the contribution to the 
power spectra from Galactic cirrus (see § 13. 3p . Each of 
the data sets is described in detail below. 

2.1. HerMES /SPIRE 

HerMES fields are organized, according to area and 
depth, into levels 1 through 7, with level 1 maps being 
the smallest and deepest (~ 310arcmin 2 ), and level 7 
maps the widest and shallowest (~ 270 deg 2 ; lQliver et akl 
12012 ). 

This study focuses on a subset of the level 5 and 6 
fields, totaling ~ 70 deg 2 , chosen for their large area and 
uniformity, and because they have a manageable level of 
Galactic cirrus contamination. The fields used for this 
study, and a summary of their properties, are given in 
Table [B The largest of the HerMES fields, the HerMES 
Large-Mode Survey (HeLMS) — which was designed 
specifically to measure the power spectrum on large an- 
gular scales — is still in preparation, and will be the 
subject of a future study. Maps will be made availa ble 
to the public through HeDaM 2 (jRoehllv et al.ll201lD as 
a part of data release 2 (DR2). 

The data obtained from the Herschel Science Archive 
were processed with a combination of standard ESA soft- 
ware and a customized software package SMAP. The 
maps themselves were then made using an up dated ver- 
sion of SMAP/SHIM (jLevenson et al.ll201Q[ ). an itera- 
tive map-maker designed to optimally separate l arge- 
scale noise from signal. SMAP differs from HIPE (lOttl 
|2010| ) in three fundamental ways which are relevant for 
power spectrum studies. First, the standard scan-by- 
scan temperature drift correction module within HIPE 
is overridden in favor of a custom correction algorithm 
which stitches together all of the time-ordered data (or 
timcstrcams) , allowing us to fit to and remove a much 
longer noise mode. Further, the standard processing is 
modified such that a "sigma-kappa" deglitcher is used in- 
stead of a wavelet deglitcher, to improve performance in 

1 http://hermes.sussex.ac.uk 

2 http://hedam.oamp.fr/HerMES/ 



large blank fields. Lastly, imperfections from thermistor 
jumps, the "cooler burp" effect, and residual glitches, 
are removed manually before map construction. De- 
ta iled descriptions of the updates to the SMAP pipeline 
of ILevcnson ct all (|2010l) ar e prese nted in Appendix |A"1 

Following lAmblard et al.l (|201lh . we make maps with 
10 iterations, fewer than SMAP's default of 20, in or- 
der to minimize the time needed to measure the transfer 
functions (§ 13.1. 2|) and uncertainties (§ 13. 2|) with Monte- 
Carlo simulations. 

Additionally, timestream data are divided into two 
halves and unique "jack-knife" map-pairs are made. 
These map-pairs are those ultimately used for estimat- 
ing power spectra (§ EH} ■ There are six sets of map-pairs 
for each field; the extra sets having custom pixel sizes so 
that the cross- frequency power spectra can be performed 
at the pixel resolution native to the maps with the larger 
instrumental beams. In other words, the maps used to 
calculate the 250 x 350 urn spectra have identical 8.33" 
pixels; and those used for 250 x 500 and 350 x 500 urn 
have 12" pixels. 

2.2. IRAS/IRIS 

At 100 urn, we use th e Improved Reprocessing of 
the I RAS Survey 3 (IRIS; iMiville-Deschenes fc Lagachd 
I2005D . a data set which corrects the original plates for 
calibration, zero level and striping problems. The result- 
ing full width at half maximum (FWHM) resolution and 
noise level are 4.3 ± 0.2', and 0.06 ± 0.02 MJy sr" 1 , re- 
spectively, and the gain uncertainty is 13.5%. Data are 
available for up to three independent observations, or 
HCONs, although two of our fields were only observed 
twice. For fields in which tiles intersect, maps can be 
stitched together with custom software provided on their 
site. 4 

3. POWER SPECTRUM OF CIB ANISOTROPIES 

The cosmic infrared background at submillimeter 
wavelengths is dominated by emission from dusty star- 
forming galaxies, while other potential sources of signal, 
like the cosmic microwave background (CMB), Sunyaev- 
Zel'dovich (SZ) effect, GHz-peaking radio galaxies or 
quasars, and intergalactic dust, are subdominant and can 
be safely ignored. Anisotropies arise from correlations 
between galaxies (i.e., galaxy clustering) which appear 
as background fluctuations. These anisotropies can be 
described by their power spectrum, and are made up of 
the following contributions: 



pshot 



pc\\ist 



of ore 
r k„ 



(1) 



Here Pj^ ot is Poisson (or shot) noise, P^ l e ust is the power 
resulting from the clustering of galaxies; i.e., the excess 
above Poisson, P[° 1C is the noise from foregrounds, and 
^ymst j g i ns + rumen tal noise. The foreground noise 
term could in principle include Galactic cirrus, free-free, 
synchrotron, and zodiacal emission, but in practice all 
but t he cirrus term are n egligible at SPIRE's wavelengths 
(e.g., lHaiian et aT]|2012D . 

The Poisson noise component arises from the discrete 
sampling of the background, and as such is decoupled 

http: //www. cita.utoronto . ca/~mamd/IRIS/IrisDownload. html 
4 http: //www. cita.utoronto . ca/"mamd/IRIS/data/irispro . tar 
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Fig. 1. — Differential (top panel) and cumulative (bottom panel) 
contributions to the Poisson (solid lines) and clustering (dashed 
line s) power from sources of different flux densities, estimated from 
the IBcthcrmin ct al.l 1120111 ) model for illustrative purposes. Top 
panel: the Poisson curves are determined by the normalization 
of S 3 dN/dS; and the clustering curves of S 4 (dN/dS) 2 . By mul- 
tiplying by an additional power of S, the peak values represent 
where the contribution to the integral per logarithmic interval are 
maximum. The integral under the curves for S cu t < 73, 61, and 
88 mJy at 250, 350, and 500 um, respectively (dotted vertical lines) 
— which represent the 5cr source detection threshold for a map 
made with five repeats — are set equal to unity. The median val- 
ues of sub-lOOmJy local maxima are (Poisson) 22, 14, and 8mjy, 
and (CIB) 7, 5, and 3mjy, at 250, 350, and 500 (am, respectively. 
Bottom panel: cumulative contribution to the power spectra nor- 
malized to unity at the detection threshold. It is evident from this 
that the Poisson level at 250 um, and to a lesser extent at 350 um, 
are very sensitive to the masking level of resolved sources, while the 
500 um channel is fairly insensitive to source masking. The clus- 
tering signal, on the other hand, is largely insensitive to masking. 
Note that the clustering curves are estimates of linear clustering 
power, and do not include non-linear effects which may be more 
sensitive to source masking. 

from the clustering term. For sources with a distribution 
of flux densities dN / dS u the effective Poisson level is 



P; 



shot 



(2) 



While the clustered power for the same sources can be 
estimated as roughly the three dimensional power spec- 
trum of the galaxy number density field weighted by the 
square of the redshift distribution of the cumulative flux, 
{dS v /dz) 2 . 

The contribution to the Poisson noise and clustered 
power from galaxies with different flux densities are il- 
lustrated as solid and dashed lines in Figure Q] respec- 



tively. The peak contribution to the Poisson noise is 
from galaxies with S v « 22, 14, and 8mJy at 250, 350, 
and 500 |J.m; while the peak contribution to the cluster- 
ing power comes from fainter (higher-z) sources, with 
S w 7, 5, and 3mJy at 250, 350, and 500 [J.m. Shown 
as dotted vertical lines at ~ 90, 75, and 60mJy are the 
5a limits of resolved sources in maps of equivalent depth 
(jNguven et al.ll2010D . Poisson noise in the power spec- 
trum is flat (in units of Jy 2 sr _1 ), behaving as a level of 
white noise which can be reduced by masking brighter 
sources. Note that masking sources is more effective at 
reducing the Poisson level at 250 |J.m than it is at 350 or 
500 iim. 

3.1. Estimating the Power Spectra: Masking, Filtering, 
and Transfer Functions 

The intensity in a given SPIRE map, / map , can be 
written as 

/map = (T ® [/sky <8> B + N] ) W, (3) 

where / s k y is the sky signal we wish to recover, T is 
the transfer function of the map-maker, B is the in- 
strumental beam, N is the noise, and W is the win- 
dow function, which includes the masking of map edges 
and of bright sources. We use ® to represent a convolu- 
tion in real-space. The instrumental noise, N , is made 
up of white noise, which dominates on angular scales 
kg ^ 0.2arcmin _1 , and 1/f noise. 

In the auto-correlation of a map (i.e., the auto-power 
spectrum of a single map), all of the power present - 
which includes both signal and noise — is correlated, 
while in the cross-correlation of jack-knife map-pairs, in 
principle the only signal correlated between them is the 
sky signal. Since we are interested in recovering the sky 
signal, we estimate the power spectrum from the cross- 
correlation of jack-knife map-pairs (discussed in § 12. 1[) . In 
practice, some correlated noise will exist between maps, 
mostly on large scales, as 1/f. To minimize this po- 
tentiality, we use map-pairs constructed by dividing the 
timcstreams in half by time, which ensures that the maps 
are made from data taken at time intervals correspond- 
ing to very large-scales. Note, this would not be true if, 
say, the data were split into those from even and odd de- 
tectors, since the same large-scale noise would be present 
in both maps. The remaining 1/f results from serendip- 
itous alignment of large-scale noise (see Appendix B of 
lHaiian et afll2012j) . 

The one-dimensional power spectrum is the azimuthal 
average of the (nearly isotropic) two-dimensional power 
spectrum of map-pairs in fc-space. In order to recover 
the true power spectrum of the sky, the cross-spectrum of 
the map-pairs must be corrected for the transfer function, 
masking, and instrumental beam. We now describe these 
corrections in detail. 

3.1.1. Masking and the Mode- Coupling Matrix 

Before calculating the power spectra, the maps are 
multiplied with windows whose values equal unity in the 
clean parts of the maps and taper to zero with a Gaus- 
sian profile (90" FWHM) at the edges. Note that jack- 
knife maps typically do not cover identical patches of sky 
because the orientation of the telescope between observa- 
tions changes, so that the windows of the two jack-knife 
map-pairs are different. 
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TABLE 2 

Number of masked sources 



Obs ID 


300 mjy 


200 mjy 


100 mjy 


50 mjy 




18 


38 


180 


1567 


BOOTES 


4 


9 


46 


723 




1 


2 


13 


248 




9 


23 


165 


1569 


CDFS-SWIR.E 


3 


g 


30 


725 




1 


1 


6 


141 




9 


23 


104 


1073 


ELAIS-Sl 


1 


3 


22 


437 




1 


1 


3 


146 




16 


42 


221 


2043 


LOCKMAN-SWIRE 


4 


9 


57 


930 




1 


3 


10 


220 




23 


66 


368 


2973 


XMM-LSS 


5 


13 


85 


1170 




1 


2 


16 


437 



Note. — Number of sources masked at 250, 350, and 500 nm 
are given from top to bottom in each panel. Not accounted for in 
the table are the extended sources masked in each field, of which 
there are 3, 2, 2, 4, and 4, respectively. A 50 mjy cut amounts to 
masking approximately 1.4, 1.2, and 0.5% of the pixels at 250, 350, 
and 500 |jm, respectively. 

Next, because we are interested in the behavior of the 
power spectrum with the masking level of resolved point 
sources, we create a series of masks for each field and 
band to mask sources whose flux densities are greater 
than 50, 100, 200, and 300 mjy, in addition to extended 
sources. Extended sources are exceptional objects, e.g., 
local IRAS galaxies, which exceed 400 mjy but can be 
as bright as 1,500 mjy at 250 urn. Mask positions are 
determined from catalogs of resolved sources which we 
construct. The source identification code first high-pass 
filters the maps in Fourier space to remove Galactic cir- 
rus and other large-scale power, then convolves the maps 
by the instrumental beam, and finally measures all peaks 
with signal-to-noise greater than 3a. The total number 
of masked sources in each field and band are tabulated in 
Table [5] Note that this me thod could potentia lly suffer 
from Eddington bias (e.g., IChapin ct al.|[2009t ). a phe- 
nomenon where instrumental white noise systematically 
boosts faint sources above the detection limit. To en- 
sure that this is not a significant problem, we compare 
the total number of s ources above t he cu t to cumulative 
number counts from iGlenn et all (|2010f l and find that 
they are consistent. 

Masks are 1.1 x FWHM in size, or 19.1, 27.7, and 
40.3" in diameter — chosen to cover the full first lobe 
of the beam, though we check that the exact size of 
the mask has a negligible effect on the spectra. For the 
data that concerns this paper, unique masks are made 
for each band, so that only sources above the given cut 
in that band are masked. This means that when calcu- 
lating the cross-frequency power spectra, not all of the 
sources masked at e.g., 250 urn will also be masked at 
350 urn, and vice-versa. However, we additionally calcu- 
late an alternative set of spectra where we mask in all 
bands the sources identified at 250 urn, meaning that the 
same sources are masked in all maps, and the power spec- 
trum pipeline is rerun. Plots and tables for this alternate 
masking scheme are presented in Appendix [Dj where we 
also show that the spectra at longer wavelengths are less 
sensitive to the level of source masking than at shorter 
wavelengths. Finally, we note that this method differs 
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Fig. 2. — Transfer functions of the map-maker, averaged in 2D, 
for each of our fields. For any given field, the transfer functions 
on large scales are indistinguishable between bands, while on small 
scales, the transfer functions converge to that of the pixel windows, 
which are band dependent. The shape of the transfer function 
depends largely on: the scan speed, which determines on what 
scales the 1// noise is projected onto the timestreams; and scan 
lengths, which determines the order of polynomial which is removed 
from the timestreams in the SMAP pipeline. Thus, maps which arc 
smaller or which were observed with slower scans are attenuated 
on smaller angular scales. 

from that of lAmblard ct all (|201lf ). who instead masked 
all pixels above 50 mjy, as well as all neighboring pixels. 
The total number of masked sources in each field and 
band are tabulated in Table [3] 

Masking in map-space can result in mode-coupling in 
Fourier space, which can bias the power sp ectrum. The 
coup ling kernel, or mode-coupling matrix dHivon et al.l 
I2002T) . in the flat sky approximation is 

Wkk'\ /N(9 k ), (4) 

where |wfcfc<| 2 is the auto-power spectrum of the mask, 
and N{9k) is the number of modes in annulus of radius 
k. Since in our case the masks of the map-pairs are not 
identical, Equation [4] is generalized for different masks 
by replacing |wfcfc<| 2 with the cross-power spectrum of 

the masks, {Wkk' w kk') (see iTristram et al"1l2005D . 

The mode-coupling matrix must be inverted in the final 
step in order to recover the de-coupled power spectrum. 
A unique mode-coupling matrix is calculated for each 
auto- and cross-frequency power spectrum, and for each 
flux cut, per field. Additionally, each Mkk' is tested on 
1000 simulated maps with steep input spectra and found 
to be unbiased. 

3.1.2. Transfer Function 

The lar ge-scale correlated n oise of SPIRE is extremely 
low (e.g.. iPascale et aT1l201ir ). As a result, only minimal 
filtering is required to make well-behaved maps, and the 
resulting power spectrum can be measured out to rela- 
tively large scales (kg ^ 0.01 arcmin -1 ). This filtering is 

quantified by the transfer function of the map-maker, T, 
which must be accounted for in the final spectra. 
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Fig. 3. — Noise levels calculated from the power spectrum of the 
difference map of jack-knife map-pairs. Not shown but consistent 
with these are noise curves estimated as the auto- minus cross- 
power spectra of the maps. White noise dominates the spectra on 
scales k$ > 0.25 arcmin" 1 , while 1/f noise is prominent on larger 

angular scales. The turnover on the largest scales, which is map- 
dependent, reflects the attenuation by the map-maker. Note the 
bump on scales kg ~ 2.2, 1.6, and 1.1, in the first, second and third 
columns, which are likely the result of their respective map pixel 
sizes (spec tra i n each column calculated from maps with same pixel 
size: see § 12.11 . 



We measure T with a Monte- Carlo simulation whose 
steps are: i) running the SMAP map-making pipeline on 
simulated observations of input maps with known power 
spectra resembling that of clustered DSFGs; ii) calcu- 
lating the power spectrum of the output map with a 
pipeline identical to that used for the real data, including 
all masking, Fourier space filtering, and mode-coupling 
corrections; and iii) computing the average of the ratio 
of the output power spectrum to the known input spec- 
trum. Transfer functions are calculated from 100 simu- 
lations of each field and wavelength, and are shown in 
Figure [2j We check that the spectra in step (i) are not 
sensitive to the steepness of the input spectra, and that 
the transfer functions have converged, with a mean error 
of ~ 1%. 

As anticipated, all three bands converge on large- 
scales for each field, meaning that the same filtering is 
performed in each band. And on small scales (kg > 
0.3 arcmin" 1 ) all fields converge to the same three curves, 
which are the pixel window functions of the three bands. 
The angular scales on which the attenuation occurs in 
each field is related to the average speeds and lengths of 
the maps scans: the former determines the scale in which 
1/f noise is projected onto the timestreams; while the 
latter dictates the order of the polynomial removed from 
the timestreams. As described later in § 14.21 the final 
spectra are weighted combinations of those in each field 
which are attenuated by less than 50%, corresponding to 
kg > 0.021, 0.009, 0.023, 0.015, and 0.009 arcmin^ 1 , for 

BOOTES, CDFS-SWIRE, ELAIS-Sl, LOCKMAN-SWIRE, and 

XMM-LSS, respectively. Maps scanned with longer and 
faster scans have less attenuation on large scales, which 



is why the maps observed in fast-scan mode (CDFS-SWIRE 
and lockman-SWIRe) are also those which best measure 
the largest scales. Note, the excess powe r introduced on 
large scales by earlier versi ons of SMAP (jLevcnson ct alj 
12010 : lAmbiard et al.ll201 l|) is no longer present. 

3.1.3. Instrument Beam 

The instrumental PSF (or beam) attenuates power on 
scales smaller than ~ 0.25-0.5 arcmin -1 , depending on 
the band. This attenuation can be corrected by dividing 
the power spectrum of the map by the power spectrum 
of the beam. The instrumental beam is measured from 
maps of Neptune — a source which to SPIRE is effec- 
tively point-like (angular size 2.5"). 

The beam power spectra are estimated in the follow- 
ing way. All pixels beyond a lOxFWHM radius from the 
peak are masked, due to uneven coverage and excessive 
noise. We check that the dependence of the beam spectra 
on the choice of the radius is small, with any differences 
contributing to the systematic uncertainties. Further- 
more, point sources in the background above 30mJy at 
250 |J.m, which are subdominant but contribute to the 
noise, are masked in all bands. Again, we check that the 
level of this masking makes a negligible difference. Thus, 
uncertainties in the beam power spectra measurements 
are largely systematic. These uncertainties couple to the 
uncertainties in the estimate of the power spectra, and 
are accounted for in the Monte Carlo procedure described 
in Sl3~2l 



3.2. Estimating Uncertainties 

Present in each map-pair is correlated signal from the 
sky, and both correlated and uncorrelated noise. Three 
terms contribute to the uncertainties in the power spec- 
trum: a non-Gaussian term due to the Poisson dis- 
tributed compact sources; sample variance in the signal 
due to limited sky coverage; and the noise. In this order, 
the variance of the cross-spectra of maps A x B can be 
written as 



° 2 (P b AxB ) = J r L + - 



p, 



AxB 



2 

, +" 

/sky n-b 

P b (N b A + N b B )+N b A N b B 



(5) 



where P b is the mean cross-spectrum of map-pairs, Nb 
is the average noise power spectrum of the map, n b is 
the number of Fourier modes measured in bin 6, and 
/sky is the observed area divided by the solid angle of 
the full sky. The first term, a\ is given by the non- 
Gaussi an part of the four-point function (as des cribed 
in e.g., lAcauaviva et al.l l2008t iHaiian et alJl2012T ). and 
is particularly sensitive to the flux cut of the masked 
sources. 

Shown in Figure |3] are the noise levels calculated from 
the power spectrum of the difference map of jack-knife 
map-pairs, which are consistent with the difference be- 
tween the auto- and cross-power spectra of the maps (not 
shown) . The noise behavior demonstrates the impressive 
performance and stability of the SPIRE instrument, with 
white noise in most cases nearly two orders of magnitude 
below the power from the sky signal. The noise spectra 
turn over on large scales due to the filtering performed 
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TABLE 3 

Galactic cirrus properties in each field 



Obs ID 


Otc 


■Po.ioo 
(Jy 2 /sr) 


-Po,250 
(Jy 2 /sr) 


-Po,350 
(Jy 2 /sr) 


B 0,500 
(Jy 2 /sr) 


T 

(K) 


BOOTES 

ELAIS-Sl 

LOCKMAN-SWIRE 

XMM-LSS 

CDFS-SWIRE 


-3.52 ±0.41 
-3.75 ±0.08 
-3.66 ±0.05 
-2.96 ±0.17 
-3.93 ±0.06 


7.22 x 10 b 
3.10 x 10 5 
2.8 x 10 5 
7.31 x 10 5 
7.36 x 10 5 


4.5 x 10 5 
1.1 x 10 6 
6.66 x 10 5 


2.23 x 10 5 

5.24 x 10 5 
3.22 x 10 5 


9.22 x 10 4 
2.9 x 10 5 
1.22 x 10 5 


17.2 ± 1.2 
18.5 ±0.7 
20.4 ± 1.4 



Note. — Best-fit variables from Equation [H] are Column 2: the index a c ; Column 3—6: the amplitudes in each band Pq; and Column 
7: the temperature T with (3 = 1.8. Long dashes represent fits which were unconstrained by data because filtering in those maps was too 
aggressive to recover the large scales where the power from cirrus would be present. 



by the map-maker, which is related to the length and 
speed of the scans. 

Uncertainties in the estimates of the power spectra are 
derived from Monte Carlo simulations of the pipeline 
on realistically simulated sky-maps. The maps include 
sources correlated between bands (i.e., the same sources 
appear in all three maps, but with different flux densi- 
ties), which are necessary for estimating uncertainties in 
the cross-frequency power spectra, as well as both 1// 
and white noise, and Galactic cirrus. Also included in 
the Monte Carlo simulations are systematic uncertainties 
arising from the beam and transfer function corrections. 

The ensemble of estimated output power spectra are 
used to measure, V, the covariance matrix 



V 



((a -a) (a-a)) mc . 



(G) 



where the tilde denotes the mean over every iteration in 
bin b. The resulting errors are 



(7 pmap 



(7) 



The non-Gaussian term emerges from the simulations as 
an offset in V. We check that this level is realistic by 
comparing it to the level of the four-point function es- 
timated directly from data, with appropriate masking, 
following IFowler et al.l (|2010D , which we found to be be- 
tween 5 and 10% of the total error in the Poisson domi- 
nated regime, depending on flux cut of masked sources: 
more aggressive source masking results in a smaller non- 
Gaussian term. 

In addition, there are ~ 8% systematic errors due to 
absolute calibration uncertainty, of which £ 1% is due 

to beam area uncertainty, as described in Appendix [Bj 
Though they are accounted for when model fitting, they 
are not included in the reported error bars. 

3.3. Galactic Cirrus 

The most significant foreground for the extragalac- 
tic power spectrum is that from Galactic cirrus, which 
can dominate the s ignal on scales greater than ~ 30'. 
iGautier et al.l (|1992h showed that the power spectrum of 
Galactic cirrus can be well approximated by a power law 



peir 



= n, I - 

fc 



(8) 



whose amplitude, Pq, normalized at fco = 0.01 arcmin -1 , 
may vary from field to field, but whose index a c ~ —3.0. 
More recently, studies of various fields using BLAST and 



e.g. q c 



—2.4 dBracco et al. 


I2011D. -2.6 dRov et al. 


dMartin et al.ll201C 


; iMiville-Dcschenes et al. 



I2010D . -2.9 ([Lagache et al.l 12007ft . and even rang- 
i ng from —2.5 to —3-6, d epending on the field 
(|Mivilfe-Deschenes et~aTI I2007D . Th us, rather than as - 
sume a power law with a = -3 fe.g.. lViero et al.ll2009D . 
we proceed by treating each field independently. 

In diffuse cirrus regions (i.e., column density, Nm < 
2 x 10 20 cm~ 2 and brightness temperature, Tb < 12 K ; 
e.g., lLockman fc Condon! [20051: iGillmon fc Slmil [2001 . 
Hi is a good tracer of dust, and the dust-to-gas ratio 
can be measured from the slope of the pixel-pixel scat- 
ter plot. Cirrus contamination can then be "cleaned" 
by scaling the Hi maps by the dust-to-gas ratio and 
subtracting them directly from the dust maps in ques- 
tion. TTris_aj)pj2acliJms_beerijK su ccessfully by 
e.g.. iPlanck Collaboration et al.l (l2011bf ) or iPcnin et al.l 
(|2012bl ). on maps which have high fidelity on large angu- 
lar scales. Unfortunately, this technique is made compli- 
cated for SMAP-made SPIRE maps because the filtering 
of scales larger than <L 20 arcmin attenuates the very 

structure that the differencing with Hi is meant to re- 
move. Though Hi maps can be filtered with the SMAP 
simulator to attenuate large scales, the remaining struc- 
tures are faint with respect to the noise, and thus difficult 
to regress with SPIRE maps. 

We i nstead adopt an approach simi l ar to that 
used by lLagachc et al.l (|2007f ). IVicro et al.l (|2009l ). and 
lAmblard et al.l (|201ll h with some additional modifica- 
tions. Diffuse Galactic cirrus emits as a modified black- 
body proportional to B{u), where B{v) is the Planck 
function and (3 is the emissivity index. Ty pically, it has 
a tem perature of ~ 18 K and (3 ~ 1.8 (e.g.. iBracco et aT] 
|2011| ) resulting in a spectral en ergy distribution (S EP) 
which peaks at ~ 170 p.m (e.g.. iMartin et al.ll2012h . At 
100 Lim, cirrus emission has roughly the same amplitude 
as at 250 urn, but unlike in SMAP/SPIRE maps, the fa- 
vorable large-scale properties of the IRIS maps make it 
possible to accurately measure the power spectra out to 
scales of ~ 4°. Thus, assuming that the Galactic cir- 
rus power spect rum is well described by a power law 
(iRov et al.M2010h . we use the 100 urn power spectra, cal- 
culated from IRIS maps with sizes identical to their 
SPIRE counterparts, and with sources above 500 mJy 
masked, to estimate the best-fit to the cirrus spectra in 
each field. Note that although larger size regions would 
better constrain the large-scale spectra, we intentionally 
use the exact same regions because these fields were cho- 
sen specifically because they were special places in the 
sky with low Galactic cirrus, and thus the spectra in- 
side and the spectra surrounding the field arc unlikely to 
be the same. Uncertainties i n the power spectra are esti- 
mated analytically following IFowler et al.1 ([20101 ). To dis- 
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Fig. 4. — Auto- and cross-frequency power spectra of the sky (circles with error bars) before correcting for Galactic cirrus (dashed lines) 
for all fields after masking sources with flux densities greater than 100 mjy. Data points are shifted horizontally for clarity. Error bars 
are derived for each field from Monte Carlo simulations, as described in § 13.21 All fields agree within the errors on scales where cirrus 
contamination is not significant ( > 0.06 arcmin" 1 ), and the variance in the shot noise levels is consistent with expectations derived from 

simulations. Maps which are smaller or which were observed with slower scans are attenuated on smaller angular scales, which is a reflection 
of the attenuation also see in the transfer functions (Figure [2} . 



tinguish between the power originating from cirrus and 
that from clustered galaxies, we include an estimate of 
the linear power (i.e., 2-hal o) term con s trained by the 
measured gal axy spectra of iPcnin et all (|20f2bf ). Also, 



we adopt the iBcthcrmin et al. (|201l[ ) model to fix the 



Poisson level, which is unconstrained by data because of 
the 4' IRAS beam, although we note that on these an- 
gular scales the exact choice for the Poisson level has a 
negligible effect on the fit. 

Next, assuming that the linear power spectrum from 
clustered DSFGs is independent of field, and after mask- 
ing all sources above 300 mJy in the SPIRE bands, we 
estimate the contribution to the SPIRE spectra from cir- 
rus by fitting the 100 \im and SPIRE auto-power spectra 
of all five fields simultaneously with: a Poisson term; a 
1- and 2-halo clustered galaxy terms; and a temperature 
(with fixed f3 = 1.8) which sets the band-to-band ampli- 
tudes, Pq- 

Lastly, uncertainties are estimated with a Monte Carlo 
simulation where the slope and amplitudes of the best- 
fit power law at 100 |xm are perturbed by an amount 
dictated by their errors, and the cirrus estimate pipeline 
described above is rerun with those values fixed. Results 
are given in Table|31 Note, gain uncertainties in the IRAS 
maps are not accounted for in the fit, as they would only 



act to increase the error of the best-fit temperature, but 
not the uncertainty in the best-fit power law. Also note 
that many of the indices are steeper than those of most 
previous analyses of cirrus power spectra, which may be 
due to the fact that these regions are specifically cho- 
sen as windows through the cirrus, and not representa- 
tive of the mean. And although they are steep, they re- 
main consistent with re sults at shorter wavelengths from 
iBazell &: Desertl (Il988f). or the extreme end of spectra 
found by iMiville-Deschenes et al.l (120071 ). Nevertheless, 
we check that fixing the slope of spectra in each field to 
a = — 3 does not significantly alter the correction, and 
indeed find that the resulting cirrus-corrected data fall 
within the uncertainties. 
We find that the method provides good constraints for 

the LOCKMAN-SWIRE, CDFS-SWIRE and XMM-LSS fields, 

while in the bootes and elais-sI fields, because of the 
aggressive filtering in the SPIRE maps, the measured 
spectra only probe scales in which the cirrus contribution 
is negligible. Consequently, when later combining the 
spectra, the largest scale bins are constrained using a 
subset of the maps. 



3.4. Re-binning and Combining Spectra 
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Individual spectra in each field, for each band and flux 
cut, are first re-binned following lAmblard ct al.1 (|2011l ) 
into logarithmic intervals with width equal to Akg /kg = 
0.25 for k e > 0.033 (A£/£ = 720), and linearly for larger 
scales, with bin-widths of Akg = 7.41 x 10~ 3 arcmin -1 
(A£ = 160). These bin-width values are chosen to ensure 
that, with the exception of correlations introduced on 
small scales from Poisson errors, the off-diagonals in the 
covariance matrix are always less than ~ 10%. The re- 
binned uncertainties are given by 



1 



(9) 



where i,j span the entries of bin b and (V^ ~) is the in- 
verse of the subset of the covariance matrix V calculated 
from simulations (§ I3.2|) . 

Next, the best estimates of the cirrus power spectra for 
each field (estimated in § 13.31 and reported in Table [3]) 
are subtracted in order to recover the power spectra of 
extragalactic sources, . Uncertainties in the cir- 

rus estimate are propagated into the final uncertainties 
assuming that the errors are uncorrelated so that 



2 2 

C pcxgal = <J psky 



(10) 



Finally, data in each band and flux cut 
are combined for the five fi elds, /, following 
iPlanck Collaboration et all ((201 lbT ) . where 



^combined 



5 

E 

/=i 



(11) 



and W[ is the weight of each field and bin, 



TP, 



/ _ 



5 

E a ps-^ 
f=i " 



(12) 



which assumes that fields are far enough apart to be un- 
correlated. Note that for each field, spectra at angu- 
lar scales where the transfer function falls below 0.5 are 
omitted in the combined fit. 

4. RESULTS 

4.1. Total Sky Spectra 

We measure signals in excess of Poisson noise in all 
auto- and cross-frequency power spectra, in each field. 
This excess signal originates from the clustering of DS- 
FGs, and to varying degrees from Galactic cirrus on large 
scales. Total power spectra of the five fields (which in- 
cludes power from Galactic cirrus and Poisson noise) 
with sources > 100 mJy masked are shown in Figure 01 
Spectra with different levels of source masking behave 
similarly, and are thus here omitted for clarity. Also 
shown are dashed lines representing the best estimate of 
the Galactic cirrus, approximated as power laws. The 
spectra from each field agree within errors on angular 
scales where the power from Galactic cirrus is subdomi- 
nant ( ^ 0.06 arcmin -1 ). 
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Fig. 5. — Spectra combined from the different fields for different 
levels of resolved source masking, plotted as circles wit h err or bars. 
The best estimates of the cirrus spectra in each field (§ [373j are re- 
moved before combining. All data and uncertainties are tabulated 
in Table [7] 
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Fig. 6. — Poisson noise level vs. flux cut of masked sources. 
Best-estimates of the auto-frequency spectra values are shown as 
open circles and tabulated in Table [8] It is worth noting that the 
data points between bands are not independent. Esti mates of the 
Poiss on level derived from SPIRE P(D) source counts I IGlcnn ct al.l 
|201"CT) are shown as exes with asymmetric error bars, whose sizes 
are functions of the uncertain upper limits on the faint end, and 
that the fields studied were relatively small. Dashed lines and 
sha ded regions repre s ent th e best estimate and lrr uncertainties of 
the lBethcrmin ct al.l d201lT ) model, with which we find good agree- 
ment at all but 350 ixm, which is undcrprcdicted by approximately 
lcr. Also shown a re measur ed Poisson levels fro m BLAST (squares: 
I Viero ct aD 12001 crosses: JHaiian et aHlWl^ ]. and SPIRE (dia- 
monds: lAmblard et aTll201lT K 
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Fig. 7. — Comparison of our cirrus subtracted, combined data to published measurements. Data are plotted as kgP(kg) in order to 
reduce the dynamic range plotted and thus better visualize the differences between the measurements. Furthermore, in order to adequately 
compare to the wide range of source masking found in the literature, we present the two masking extremes of our analysis: spectra with 
sources greater than 50mjy maske d (dark blue cir c les), and those with only extended sources having been masked (light blue circles). 
Previous SPIRE measurements from lAmblard ct al.l l|2011l ). which should be compared to the dark blue circles, are shown as black crosses. 
The remaining data and curves should be compared to the light blue circles. Note t hat to help a id the comparison, error bars include 
system atics due to calibration and beam uncertaint ies. They are: BLAST data from l Vicro ct ahl I I2009I ) (brown exes) and IHaiian ct aTl 
(120121 ) (lavender diamonds); and Planck data from [Planck Collaboration ct ahl J2011bl ) (red squares) at 350 and 500 Ltm. Note, Planck 
data at 857 and 545 GHz (350 and 550 Ltm) arc color corrected to account for their different passbands by multiplying them by factors of 
0.99, and 1.30, respectively. 

iBethermin et alj (1201 ll ) . a phenomenological model 
which to-date is the best at reproducing the observed 
number counts from 15 |xm to 1.1mm. We find that the 
model is in very good agreement with the data at all flux- 
cut levels and at all but 350 |xm, which underpredicts the 
data by approximately la. 



4.2. Combined Extragalactic Sky Spectra 

The cirrus-subtracted power spectra in each field, 
which are combined using Equation II H are presented in 
Figure [5] and tabulated in Appendix [Ej Table [71 As ex- 
pected, the spectra in each panel converge for increasing 
angular scales as the contribution from Poisson noise be- 
comes subdominant. Consistent with expectations from 
Figure [TJ the flux cut has a significant effect on the Pois- 
son level at 250 urn, and a nearly negligible effect at 
500 Ltm. 

Poisson levels are determined through a simultaneous 
fit to the combined spectra of the Poisson and clus- 
tered galaxy term s with templates adopted from the 
IViero et al.l (|2009l ) halo model, and are shown as a func- 
tion of flux cut of masked sources in Figure [5] We 
note that this estimate is subject to systematic un- 
certainties due to the mild degeneracy of the Poisson 
and 1-halo terms, and that we account for those uncer- 
tainties in the estimate. As anticipated, shorter wave- 
lengths are significantly more affected by removal of 
the brightes t sources (see F i gure [p . Previous measure- 
ments from lAmblard et al.l (|2011l ) with 50mJy sources 
masked are in relatively good agreement, with 250 urn 
higher by ~ 5%, and 350 and 500 Ltm lower by ~ 7% 
and 17%, respectively. The BLAST measurements, 
which m asked sources with flux densities gr eater than 
500 mJy (|Viero et al.ll2009HHaiian ct al.ll2012l ) appear to 
be higher than our values by ~ 26%, 12%, and 11% at 
250, 350, and 500 Ltm, respectively. 

Also pl otted are est i mates of the Poisson level derived 
from the iGlenn ct al.1 (l2010f ) P(D) number counts using 
Equation [2] We find that the number count predictions 
over-estimate our values by ~ 16% at 250 |xm, and under- 
estimate our measured values by ~ 12 and 16% at 350 
and 500 |J.m, respectively. 

Finally, we compare to the model predictions of 



4.3. Comparison to Published Measurements 

In Figure[7]we plot our combined auto-frequency power 
spectra along with a selection of recently published CIB A 
measurements. We show the two masking extremes of 
our data: those in which all sources greater than 50 mJy 
were masked (dark blue open circles); and those where 
only extended sources were masked (light blue open cir- 
cles). We do this in order to adequately compare with 
the wide range of masking in the literature. 

Shown as gray points are the SPIRE auto-frcquenc 
power spectra of ~ 15deg 2 from lAmblard et al.l (|201ll 
in which pixels greater than 50 mJy were masked, so that 
they should be compared to our dark blue circles. On 
larger scales (e.g., kg ^ 0.08 arcmin -1 ), particularly at 

shorter wavelengths, their spectra suffer from overcor- 
re ction of Galactic cirrus contamina tion (see discussion 
in lPlanck Collaboration et al.ll20lTbt ). On smaller scales 
we find that our spectra differ by factors of ~ 0.91 ±0.01, 
1.09 ± 0.01, and 1.09 ± 0.01. Note that the calibration 
of the science demons tration phase (SDP) maps used in 
lAmblard et "all (pOll differed by 1.02, 1.05, and 0.94, 
at 250, 350, and 500 |J.m, respectively, but that those 
corrections were not applied here. These calibration dif- 
ferences, combined with the offsets resulting from the 
new estimate of the beam (Appendix |Bj may partially 
account for this difference. Also note that although the 
error bars on their data are comparable to ours at small 
angular scales, correlations due to the non-Gaussian term 
(first term in Equation [SJ were not included when they 
re-binned into log bins, thus artificially deflating their 
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Fig. 8. — Comparis on of our cirrus-corrected, combined data to published models. Best-fit halo model to previously published SPIRE data 
IIAmblard et al.N20I II ) are shown as grey dot-dashed lines, and should be compared t o the dark blue c ircles. The remaining curves should 
be compared to the light blue circles. Note, with the exception of the BLAST model llVicro ct al,H200"9T) , all of the following were originally 
fit to the 857 and 545 GHz chann els of Planck, and have thus been color-corrected by .99 and 1.30 at 350 and 500 |-im, respectively. 
They are from: IVicro ct al.l (120091 . brown do tted lines) ; [P lanck Collaboration ct all (|2011bl . red dashed lines) at 350 and 500 nm; case of 



IShang ct all (I2012I . green dashed lines); and IXia ct all (120121 . orange three-dot-dashed lines) 



errors. 

Shown as red squares at 350 and 500 \wi are re- 
sults from the lPlanck Collaboration et alj (|2011bD . which 
should be compared to our light blue circles. Note, com- 
parisons of the two spectra must be made with caution, 
bearing in mind that the flux density of masked sources 
in Planck is much higher (710 and 540 mJy at 350 and 
550 urn) . In addition, because the passbands of the two 
instruments are not the same, Planck data at 857 and 
545 GHz (350 and 550 |a.m) are color corrected by multi- 
plying them by factors of 0.99, and 1.30, respectively. 

We find possible inconsistencies between the two sets of 
measurements. On scales kg ;$ 0.04 arcmin" 1 , the Planck 

spectra appear to be offset high by factors of 1.26 ± 0.06 
and 1.40 ± 0.06 at 350 and 500 u.m, respectively; while 
on smaller angular scales there is an apparent excess of 
power in the Planck data, though it should be noted that 
point-by-point, the 350 |J.m values do agree within errors. 
Potential explanations for this discrepancy include cali- 
bration, excess Poisson, or reconstruction systematic er- 
rors in the Planck beam. We discuss these scenarios in 
more detail in § 16.21 



4.4. Comparison to Published Models 

In Figure[5]we plot our combined auto- frequency power 
spectra next to a selection of published halo models. As 
in Figure [3 we show the two masking extremes of our 
data in order to adequately compare with the wide range 
of mas king in the l i teratu re. 

The IViero et al.l (|2009D models (dotted br own lines), 
which were fit to BLAST data with sources greater than 
500 mJy masked, and appeared to be a good match to 
the iPlanck Collaboration ct all (|2011bl ) data, here do a 
poor job of describing the SPIRE measurements, overes- 
timati ng the power on scales greater than ~ 40". 

The lAmblard et alj (|201lD models, which assumed a 
masking level of 50 mJy and were fit to each band indi- 
vidually, are shown as grey dot-dashed lines. Similar to 
the differences in the data, the large scale power is under- 
estimated due to the over-correction for the contribution 



from Galactic cirrus; while on small scales, the small dif- 
ferences may be due to diffe rence in the calibration. 

The h alo models from the IPlanck Collaboration et al] 
(l2011bl ) are shown as red dashed lines at 350 and 500 |xm. 
They assumed a masking level consistent with the mask- 
ing level of Planck data, as do all following models. They 
agree very well with our data on kg ^ 0.04, at which 

point they turn more sharply than the data, leading to 
excess power on smaller angler scales. We discuss possi- 
ble rea sons for th i s diffe rence in § 16.21 

The IXia et al.l (|2012l ) halo model, shown as orange 
three-dot-dashed lines in Fig ure El was fit to Planck and 
correct ed SPIRE data from IPlanck Collaboration et ail 
(|2011bl § 5.3). It adopts a descr iption for the 
source population fr om lLapi et all (|201ll ) (an update of 
iGranato et al.l I2004I ). It appears to be consistent with 
the overall amplitude of the data, but has a bump of 
excess emission at around kg ~ 0.1-0.03 that the data 
do not show. This evolutionary model assumes that the 
steep part of the source counts at submillimeter wave- 
lengths is dominated by massive, proto-spheroidal galax- 
ies in the process of forming most of their stars. The 
model also includes small contributions from late-type 
and starburst galaxies. Notable in this model is that 
the redshift distribution of the emission peaks at slightly 
higher redshifts, broadly around z ~ 1.7-2.2, increas- 
ing with increasing wavelength; this is in distinction to 
other models which are strongly peake d at z ~ 1. It 
was fit to data from ffersc/teZ/SPIRE (lAmblard ct all 
120111). Planck /RFl IPlanck Collaboration ct al.l (l2011bD 



SPT (jShirokoff ct all 1201 ll ) and ACT (jDunklev ct all 
120111 ). and as it is physically based, it is much more con- 
stra ined than the phe no menological models used by of 
e.g.. lViero et~all (|2009D or lShang et all (|2012D . That two 
populations are represented is evidenced by the clear fea- 
ture at ~ 0.2 arcmin -1 , a feature which is not apparent in 
the data, suggesting an overestimate of the contribution 
of late-type galaxies to the total spectrum. 

Las tly, shown as green dashed lines is the lShang et al.1 
(|2012f ) model, whose main feature is to implement a 
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Fig. 9. — Clustering spectra combined from different fields vs. flux-cut leve l of masked sources . Spectra are shifted horizontally for visual 
clarity. Overlaid are best-fit templates to the data, which were taken from Ivlcro ct al] d2009r) . If resolved sources contributed solely to 
the Poisson noise component of the spectra then these points would lie on top of one another. Instead, there is a clear reduction of 1-halo 
power with masking level at 250 |jm, suggesting that some fraction of bright DSFGs reside, as satellites, in more massive dark matter halos. 



luminosity-mass {L — M) relation, such that more mas- 
sive halos host more luminous sources. Though the 
model was fit primarily to Planck data, it appears to 
fit our spectra at 500 u.m quite well. The fit is less good 
at 250 and 350 u.m; though the shape is in good agree- 
ment, the curves are high by ~ 30%. Despite this success, 
the model has some points of concern. In particular, it 
underpredicts the contribution of lower redshift sources 
to the CIB (e.g., they are significantly below the lower 
limits measured f rom stacking 24 u.m selected sources; 
iJauzac et al.ll20ilh , and also find an uncharacteristically 
large j3 (discussed further in § 16. 3|) . 

Finally, we highlight the feature that appears in the 
data at kg ~ 0.03-0.04 arcmin -1 , particularly in all 
auto- and cross-frequency spectra, but predominantly at 
50 u-m. A similar featur e was visible in the SPIRE data 
of lAmblard et all (|2011[) as well. It is not clear if this 
is a real feature in the sky, which would be unexpected 
and is unlikely, or due to noise present in the data. Sim- 
ilar features do not appear in the transfer function, or in 
the simulations of the pipeline which tested for potential 
biases. 

4.4.1. Clustered Galaxy Power Spectra 

Ultimately, we arc interested in the power spectra 
of clustered DSFGs, which we estimate by removing 
the Poisson noise from the cirrus-subtracted, combined 
power spectra of Figure[5] Results are shown in Figure[9l 



The fraction anisotropy, calculated as 
SI /I = 2irk 2 g P(k e )/L 



GIB, (13) 

where Icib is the overall amplitude amplitude of the 
CIB, measured to be 0.71 ± 0.17, 0.59 ± 0.14, and 
O.SSd iO.lOMJysr- 1 at 250 350, and 500 urn, re spectively 
fe.g. jLagache et al.|[2000HMarsdcn et al.ll2009ft. We find 
SI/I = 14±4%, consistent with findings from lviero et al.1 
(l200l . 

We fit a simple power law to the clustering spectra 
over the range kg = 0.01-1.4 arcmin^ 1 , finding a mild 
change in the slope with changes in the masking level. 
Specifically, from most (S^ut > 50mJy) to least aggres- 
sively masked (only extended sources), we find slopes of 
-1.60 ±0.05 to -1.50 ±0.07 at 250 u.m; -1.52 ±0.05 to 
-1.36±0.06at 350 urn; and -1.52±0.06 to -1.47±0.06 
at 500 u.m. The x 2 of these fits for 15 degrees of freedom 
(d.o.f.) are 7, 10 and 3 (reduced \ 2 ~ 0.5, 0.7, and 0.3) 
at 250, 350, and 500 u.m, respectively. 

The next notable feature is the reduction of 1-halo 
power with flux cut of masked sources, particularly at 
250 u.m, which are shown as dashed lines in Figure [9j 
whereas the 2-halo power, shown as dotted lines, remain 
relatively unchanged. We demonstrate in § 16.11 how this 
result can be interpreted as more luminous sources re- 
siding in more massive halos, motivating the use of a 
model later in the paper in which a luminosity-mass rela- 
tionship is invoked (e.g.. IShcthll2005t ISkibba et alj|2006t 
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Fig. 10. — Cross-correlation power spectra. Spectra are shifted 
horizontally for visual clarity. For two identical maps, or two maps 
in which all the sources are at the same redshift, this measurement 
would have unit amplitude, which is represented by a dotted line. 
We find that for unmasked maps, the cross-correlation is approx- 
imately 0.95 ± 0.04, 0.86 ± 0.04, and 0.95 ± 0.03, for 250 x 350, 
250 X 500, and 350 X 500, respectively. For maps with sources 
greater than 50 mJy masked, this cross-correlation is reduced, and 
appears to weaken further with decreasing angular scale. 



IShang ct al.ll2012T ) . We show that attempting to account 
for the reduction of power entirely with the Poisson term 
leads to significant tension in the fit. Though the reduc- 
tion of power is much less significant at 500 u.m, we re- 
mind the reader that there are far fewer sources at each 
given flux-cut level at 500 u\m than at 250 or 350 [J.m (see 
Table EJ. The capability to mask fainter sources reliably 
at 500 urn would require either maps with higher angular 
resolution (i.e., less confusion noise) or a way to probe 
deeper into the confus ion using ancillary data (e.g., XID; 
iRoseboom ct al.ll2010D . 

4.5. Cross- Correlation Power Spectra 
The cross-correlation power spectrum is defined as 



Caxb — 



P, 



AxB 



fcfl 



(14) 



i.e., the ratio of the cross- frequency power spectra to the 
geometric mean of the two auto- frequency power spectra. 
Identical maps would thus have a cross-correlation of unit 
amplitude, as would maps containing sources located at 
identical redshifts and with identical colors. Departures 
from unity would be an indication that sources are not 
all at the same redshift, or that their colors (or average 
temperatures) are variable, and the strength and shape 
of the cross-correlation signal would depend on the level 
of correlation between maps. Consequently, the cross- 
correlation provides strong constraints for source popu- 
lation models. 

We show the cross-correlations as functions of the 
flux cut of masked sources in Figure [K3 The measure- 
ment becomes very uncertain at angular scales kg < 
O.laxcmin . At larger kg, with the exception of the 



50mJy cut, we find similar levels of correlation for all 
levels of masking, which can be approximated as hori- 
zontal lines at 0.95 ± 0.04, 0.86 ± 0.04, and 0.95 ± 0.03 
for A x B = 250 x 350, 250 x 500, and 350 x 500, respec- 
tively. Cross-correlations of maps with sources masked 
at 50 mJy appear to be less correlated, with hints of a 
reduction in the c orrelation wi t h incr easing kg, which as 
first predicted by iKnox et all (|200lh , would be an in- 
dication that longer wavelengths are more sensitive to 
higher- z. 

These results compa re favor ably with the 
iPlanck Collaboration ct ail (|2011bO measurements 
of the cross-correlation, who found 0.89 and 0.91 for two 
different fields at 350 x 550 (i.m, and are also consistent 
within errors with the cross-correlations measured by 
lHaiian et all (|2012ft . 

5. HALO MODEL INTERPRETATION OF CIB 
ANISOTROPY MEASUREMENTS 

The angular power spectrum of intensity fluctua- 
tions , P uu >(kg), is obtained using Limber's approxima- 
tion (lLimbcrlll953f ) , which is valid on small angular scales 
(2irke /t 10). The projection of the flux- weighted spatial 
power spectrum is 



_ ., . f dz dz /, 
Puw(k e ) = J ^^ p ^' ( k 



lixkg 



dS v dS„i 
dz dz 



(15) 

where x( z ) lS the comoving distance to redshift z, and 
dSyjdz is the redshift distribution of the cumulative flux. 
For sources with flux densities S v < S cv x 



dS v 

dz 



d 2 N 
dS v dz 



{S v ,z)dS v , 



(16) 



and the differential number counts are related to the 
epoch-dependent comoving luminosity function, <f>(L,z), 
through 



dN 
dS v 



,2 d X 



dzx -r4>[L{S„,z),z] 
dz 



(17) 



5.1. Halo Model Formalism 



The power spectrum of CIB anisotropies in the halo 
model formalism is written as the sum of three terms: 
the linear (or 2-halo) term, which accounts for pairs of 
galaxies in separate halos and dominates the spectrum 
on large scales; the non-linear (or 1-halo) term, which 
describes pairs of galaxies residing in the same halo and 
is the dominant term on small scales; and the Poisson 
(or shot) noise term: 

P vv , (fc, z) = Pll, (k, z) + Pll, (fc, z) + P*? (fc, z). (18) 

Common to most halo models is that a distinction is 
made between central and satellite galaxies, with A^ gal = 
N ccn + 7V sat . All halos above a minimum mass M m - m host 
a galaxy at their center, 



N c 



v , r /0 M<M min , 
( ' M M>M min , 



(19) 



while any additional galaxies in the same halo would be 
designated as sa tellites which trace the dark matter den- 
sity profile fe.g. JZheng ct al1l2005f ). Halos host satellites 
when their mass exceeds the pivot mass Mi (also known 
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as M sat in the literature) , and the number of satellites is 
an exponential function of halo mass: 



N sat (M) 



M 
Mi 



(20) 



In the standard halo model, galaxies are as- 
sumed to c ontribute equally t o the emissivit y den- 
sity (e.g., IViero et all l2009t lAmblard et all 120111 : 
iPlanck Collaboration ct alj|2011br ). Therefore, assuming 
that the CIB originates from galaxies, spatial variations 
in the specific emission coefficient j v directly trace fluc- 
tuations in the galaxy number density 

Sju/ju = Sn^/n^ 1 . (21) 

The linear, 2-halo term dominates on large scales and 
is given by the clustering of galaxies in separate dark 
matter halos: 

P%>(k,z) = ^L T P Xin {k,z)D v {k,z)D' v {k,z), (22) 



3u3 



with 



D v {k,z) = J dM 



dN 
dM 

x[N^ cn (M, z) 



(z)b(M,z)u sa i(k,z\M) 
K at {M,z)}, 



(23) 



where Pn n (k, z) is the linear dark matter power spec- 
trum, 6(M, z) is the linear large-scale bias, and 
u g ai(k, z\M) is the normalized Fourier transform of the 
galaxy density distribution within a halo, which is as- 
sumed to equal the dark matter density profile, i.e., 
u ga \(k,z\M) = UY>yi{k, z\M). The non-linear, 1-halo 
term dominates on small scales, and is written as 



P^(k,z)- 



/ dM—(z) 
j v ji J dM y > 

x N^(M\z)N^(M\z)u gal (k,z\M) + 
N c ™(M\z)N: at (M\z)u gal (k, z\M) + 
N s f{M, z)N° a \M, z)u 2 gal (k, z\M)\ , (24) 



(jCoorav fc Shethll2002l) . where dN/dM is the halo mass 
function. 

However, conceptually it is wrong to assume that 
galaxies of different luminosities have equal weight in 
contributing to the power spectrum of the intensity fluc- 
tuations. A consequence of this assumption is that 
the excess signal on small angular scales from galax- 
ies in massive halos can only be reproduced by having 
more satellite galaxies, leading to previous estimates of 
a which exceed predictions for sub-h alo indices from 
semi-analytic mod els, (a < 1: e.g.. iGao et al.l [20041 : 
lHansen et al.l l2009t ) and a significant overabundance of 
satellites. For example, pre vious fits to SPIRE dat a with 
a standard halo model from lAmblard ct al.l (|201ll ) found 
a ~ 1.7 at 250 urn, and ~ 1.8 at 350 and 500 |J.m, albeit 
with large errors. 

A way to overcome this excess satellite problem, while 
still producing enough 1-halo power, is to have a model 
with fewer but more luminous satellites and weight galax- 
ies by their l umin osities. One such model is that of 
IShang et alj (|2012f ). which invokes a luminosity- mass 
(L — M) relation to tie the emissivity from galaxies to 



their host halo masse s (also see e.g., lYang et al.l 12001 
I Vale fc Ostrikcrl 12001 ). The advantage of this model is 
that in principle one can predict the abundance as well as 
the clustering of galaxies observed at different frequency 
bands simultaneously, while in the standard halo model 
it is impossible to predict the power spectrum across dif- 
ferent frequency bands at the same time. As we will now 
show, this new formalism is similar to the standard one 
up until the number of central and satellite galaxies are 
substituted for luminosity weighted quantities. 

5.2. Luminosity Weighted Halo Model 

Her eafter, we follow the formalism of IShang ct al.l 
(|2012t ). with some modifications. Novel to our implemen- 
tation is the simultaneous fitting to the power spectra in 
all three bands, and to the number co unts of sources 
above about 0.1 mJy (jGlenn ct al.l l2010h . In this imple- 
mentation of the halo model, the mean comoving specific 
emission coefficient is 



3v{z) 



dn . L v [{l + z)v] 

~d~L 4^ ' (25) 



where L is the luminosity and dn/dL is the luminos- 
ity function of DSFGs. What this means is that, unlike 
the standard model which assumes that all galaxies con- 
tribute equally to the emissivity density (i.e., have the 
same luminosity), here the emissivities of galaxies in a 
given halo, by way of their luminosities, depend on the 
redshift, halo mass and frequency: 

L„[(l + z)v] = L (l + z)"S(m)9 [(1 +z)v], (26) 

where Lo is the overall normalization factor, rj describes 
the redshift evolution, S(m) describes the relation be- 
tween infrared luminosity and halo mass (the L — M 
relation), and 0(u) describes the shape of the infrared 
SED. 

Furthermore, the effective, luminosity- weighted, num- 
ber of central and satellite galaxies is 



fr = at 



t cen 
cen (1+2) 



v (M,z) 



4tt 



J v 



dm— — (M, z) 
dm 



L 



{l+z)v 



(m, z) 



(27) 



(28) 



47T 

The terms in Equations [521 an d HH remain the same, 
except the number of central and satellite galaxies are 
substituted with their luminosity weighted counterparts: 
such that D v in the 2-halo term (Equation |2"31 becomes 

dN 



D u (k,z) = J dM 



dM 



(z)b(M,z)u ga i(k,z\M) 



x[r™{M 1 z) + fl at {M,z)] 1 
and the 1-halo term (Equation l24l) becomes 
1 f ,..dN 



(29) 



dM ^) 



x /£f n (M, z)f^(M, z)u mX {k, z\M) + 

f c r{M, z)ff{M, z)u sai (k, z\M) + 

fi?(M, z)f?(M, z)u 2 gal (k, z\M)] . (30) 

We define halos here as overdense regions whose mean 
density is A„^ r times the mean density of the Universe 
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according to the sph erical collapse model, and we ad opt satellite galaxies. Equation! 



the density profile of lNavarra Frenk. k White] (Il997l. 



NEW) with the concentration par ameter of 



or 



Bullock ct all 



(|200lft . and the fitting function of lTinker et al.l (|2008l) for 



the halo mass fu nction and its associated prescri ption 
for the halo bias (jTinker. Wechsler. k Zhend[2010l ). For 
the subhalo mass functi on, we use the fitting function of 
iTinker k Wetzell (|2010ft . We will now describe each of 
the terms in Equation [2U in more detail. N ote that us- 
ing in stead the concentration parameter of iDuffv et alj 
(|2008l ) leads to very little change in the final best-fit pa- 
rameters. 

5.2.1. (1 + z) v : The Luminosity Evolution 

The luminosity evolution in this model is motivated by 
the known increase of specific star- formation rate (sSFR) 
with redshift (e.g. lElbaz ct alj|2007t lOliver ct al1l2010at 
Karim et all 1201 It iNocskc et all l2007t ISargcnt et all 



20121 IWang et al.ll2012l ). and the fact that SFRs and 
infrared lumino sities are correlated for DSFGs (e.g., 
iKennicuttl 11998ft . The exact form of the evolution is 
still not clear: though measurement s find a rap i d rise 
followed by plateau at z £ 2 (e.g., IStark ct"all 120091 : 

iGonzalez et all I2010T ) , semi-analytic models have diffi- 
cultly reproducing observations without invoking a num- 
ber of ad-h qc modifications to the standard p hysical 
recipes (e.g. IWeinmann. Neistein. k Dekcll 1201 If ). Yet, 
without a convincing alternative, we proceed motivated 
by observations, letting 77 be a free parameter over < 
z < 2 and set 77 = at z > 2. We also explore the 
scenario where a single 77 is applied to the entire redshift 
range considered (0 < z < 4) and find that it gives a 
poorer fit to the data. 

5.2.2. E(m): The L — M relation 

Observationally it is clear that some halos are more 
efficient than others a t hosting sta r form ation (e.g., 
iBethermin et all l2012at IWang et al] 120121 ). and that 
the halo mass of most efficient star formation evolves 
with redshift (i.e.. downsizing; e.g.. ICowie et all 119961 : 
iBundv et all I2006T ) . It is also clear that star forma- 
tion in halos is suppressed by several plausible mech- 
anisms at the high mass (e.g., accretin g black holes 
iBirnboim k Dckelll2003l : iKcres ct al.ll2005ft and low mass 
(e.g., feedback from supernovae, photoiqnizat ion heating 
iDckcl k Silkl [19861: iThoul k Weinberg] fl996h extremes. 
Thus, following IShang ct al. I (120121 1. we assume that the 
L — M relation, S(m), can be parameterized by a simple 
log-normal distribution 



S(m) 



2itct 2 t i 

L 1 m 



■. exp 



(logm - logTOp Cak ) 2 

L/rn 



(31) 

where rn poa k describes the peak of the specific IR emis- 
sivity per unit mass, and 0£/ m describes the range of 
halo masses in which galaxies producing IR emission re- 
side. The minimum halo mass to host a galaxy, M m - ln , is 
left as a free parameter, but we place a lower limit on it 
such that L = at M < M min . 

Note that we have implicitly assumed that the shape of 
the relation between halo mass and infrared luminosity is 
rcdshift-independent and identical for both central and 



, , , , dN , . 1 

jjz)= / dM (z) — 

JuK ' 1 dM y Uir 



can then be recast as 

^ycen £ccn 



(171 

dm—(M,z)L sat 
dm 



{l+z)v 



(32) 



where m is the subhalo mass at the time of a ccre- 
tion (e.g. JWetzel k Whitdl201(il : IShang et al.ll2012h and 
dn/dm is the subhalo mass function in a host halo of 
mass M at a given redshift. 

5.2.3. Q(v): The Model SED 

Following lHall et all (I20T0I ) and IShang et all (pOll . 

we adopt a single modified blackbody SED, where 

@{v)(x^B{v,T A ), (33) 

where B(v,Td) is the blackbody spectrum (or Planck 
function), with effective dust temperature Td, and /3 is 
the emissivity index. Both Td and j3 are left as free vari- 
ables. Similarly, we replace the mid-infrared exponential 
on the Wien sid e of the spectrum with a power- law of the 
form /„ a ' a (|Viero et al.ll2012l ). which we fix to a = 2. 
This technically requires that the two functions and their 
first derivatives are equal at the transition frequency. 

5.2.4. Markov Chain Monte Carlo 

We make use of Markov Chain Monte Carlo (MCMC) 
methods to derive the posterior probability distributions 
for all parameters by fitting to the P(D) number counts 
dGlcnn et al.ll2010T ) and auto-frequency power spectra at 
250, 350 and 500 um simultaneously. Due to computing 
time considerations of the MCMC process, we fit to only 
the spectra with sources above 50 and 300 mJy masked. 
Furthermore, we do not fit to the cross-frequency power 
spectra, but rather we predict the cross-frequency power 
spectra using the best-fit model parameters, and then 
compare with the measured cross-frequency power spec- 
tra. Our model consists of seven parameters, consisting 
of: two for the SED (T d and /3); two for the L - M re- 
lation (Mp 0a k and c\L/ m ); one for the redshift evolution 
(77); one for the low mass halo cutoff (Af m ; n ); and an over- 
all normaliza tion (Lp ). As some of these are degenerate 
IShang ct all (see e.g. . 120121 1 we later show that Td and rj 
can reasonably be fixed, leaving five free parameters. In 
comparison, there are a total of eight pa rameters in the 
standa rd halo model as implemented in lAmblard et all 
(|201ll ). and moreover, the parameters in the standard 
halo model in each frequency band must be derived sep- 
arately, reducing the constraining power of the data. 

5.3. Halo Model Results 

We first explore a model with all seven parameters free 
and find significant correlations (of order unity) between 
several of them, particularly the pairs: Td~/3, M pB3 ^-rj, 
and Lq-t]. Thus, we perform the fit again for two more 
cases: A) first fix ing Td = 34 K and 77 = 3.25 (following 
IShang et al]|2012f l: and B) fixing /3 = 1.8 lDraine fc Led 
(e.g.. 11984ft ' and 77 = 3.25. Results of all three cases are 
tabulated in Table H The X 2 = 256 for 130 d.o.f. (or 
Xreduced = 2 -0) m tne ^ parameter fit, and \ 2 = 250, 245 
for 125 d.o.f. (or Xreduccd = 2-0, 2.5) in the 5 parameter 
fits with T and (3 fixed, respectively. Further, for case 
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Fig. 11. — Best-fit hal o model fit s i multa neously to cirrus-subtracted, combined auto-frequency power spectra (top row) and to the 
P(D) number counts from [Glenn ct aTl 1120101 ) (see Figure [T2j . Spectra where Scut = 300 mJy (upper green circles) and Scut = 50mjy 
(lower blue circles) are fit with three terms: Poisson (horizontal dashed lines); 2- halo (steep dotted lines); and 1-halo (nearly horizontal 
dotted lines). The sum of the three terms are shown as solid lines. Plotted in the bottom row panels are the best-fit points minus the data, 
divided by the data. A perfect fit would have a null value, represented here as a horizontal dotted line. 



TABLE 4 
Best-Fit Halo Model Parameters 



Case r, T [K] g log(M min ) [M ] log(Af pcak ) [Mp] aj /m log(L ) [L ] X r 2 cduccd 

7 Parameter 3.16 ±0.25 34.6 ±2.1 0.77 ±0.13 10.9 ± 0.7 13.8 ± 0.6 0.29 ± 0.02 -2.6 ±0.1 2.0 

5 Parameter A (3.25) (34) 0.83 ± 0.04 10.8 ± 0.6 13.9 ± 0.6 0.28 ± 0.02 -2.6 ±0.1 2.0 

5 Parameter B (3.25) 25.5 ± 1.1 (1.8) 10.9 ± 0.6 13.9 ± 0.6 0.24 ± 0.01 -2.6 ±0.1 2.5 

Note. — Best-fit values and marginalized la errors from the combined fits to the spectra and source counts. The first row corresponds 
to a fit where all seven parameters are allowed to vary, which led to significant degeneracies. In the next two rows, the parameters r\ and: 
T in Case A; /3 in Case B, are fixed, which are indicated by their placement in brackets. Lastly, the \ 2 in the last column correspond to 
130 degrees of freedom for the 7 parameter fit, and 125 degrees of freedom for the 5 parameter fits. 



TABLE 5 

Correlations between best-fit parameters of halo model 



Parameter 


T 




M min 


-^pcak 


° 2 r, 

L/m 


Lo 


T 

Mmin 


1 


1 

0.015 


-0.001 


-0.098 


-0.031 


0.097 




1 


-0.010 


0.009 


-0.011 


-^peak 




-0.586 


-0.009 


1 


-0.009 


0.334 


L/m 




0.080 


0.014 


-0.002 


1 


-0.767 


Lo 




0.567 


0.008 


-0.071 


-0.323 


1 



Note. — Correlations between parameters represent degenera- 
cies of the fit. The lower triangular are values for the 5 parameter 
fit with T fixed (Case A in Table 4), while the upper triangular is 
the same but with (3 fixed (Case B in Table 4). Values of +1 or — 1 
mean that the parameters are highly correlated or anti-correlated, 
respectively, while values near mean that they are independent 
of one another. 

A (fixed T), the counts and the spectra are equally well 
fit, with reduced \ 2 ~ 2.0. Hereafter, we will refer to 
the case with T fixed, as it provides a superior fit, unless 
otherwise noted. 

Our best-fit halo model is shown in Figures HT1 and [T2l 
Also shown are predictions from the model for the cross- 



frequency power spectra in Figure [T3l The spectra are 
remarkably well fit by the model, with both the Poisson 
and 1-halo terms reducing in amplitude, while the 2-halo 
term remains relatively stationary. The counts are less 
well fit, which may reveal a tension between fitting the 
spectra and the counts together, or that perhaps more 
freedom (particularly in the SEDs) is needed. 

The rcdshift distribution of the emissivity — which 
in the standard halo model is either parametrized or is 
adopted from galaxy population models — is here an out- 
put of the L — M relation, and shown c ompared to a se- 
lectio n of models (jValiante et al.ll2009t iBethermin et al.l 
1201 It Bother mi n ct a l. in prep.) and previous estimates 
(jAmblard et al.ll2011[ ) in Figure 1141 The behavior of our 
rcdshift distribution is similar to the models, in that they 
predict a fast rise from z = to 1 or 2, but disagrees at 
higher z, where our model contribution declines much 
more rapidly. It should be noted that the models at 
z 2 remain highly uncertain. 

We find that star formation is most efficient in halos 
ranging from log(M Q ) = 13.3 to 14.4, consistent with the 
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between the standard model and ours. 
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Fig. 12. — Euclid ean normalized differential number counts from 
I Glenn et al.l 1 120101 . shown as circles) at 250 (blue), 350 (green), 
and 500 \xm (red) along with best-fit curves from the model (solid 
lines). These curves were found by simultaneously fitting to spectra 
(shown in Figure [TTJ as well as to these counts. This fit reveals 
a possible tension between the modeling of the counts and the 
modeling of the clustering terms. 

halo masses for galaxies which become b ulge-dominated 
galax ies in semi-analytic models (e.g., iGonzalcz et al.l 
1201 If ). However, our Mpoak is much higher than the 
McS ~ 5 x 10 12 Mm for z ~ 2 sub-millimeter g alaxies 
found by IViero et ail (|2009j) and IXia et all (|2012f) . It is 
worth noting that in the standard halo model, the ef- 
fective halo mass scale is calculated as the average halo 
mass weighted by the number of galaxies in the halo 



M, 



elf 



dn A^ al 
dM — M——r- 

dM (iVs al ) 



(34) 



As such, M c s should not be directly compared with 
Afpcak in our model, which corresponds to the halo mass 
scale where luminosity-to-mass ratio is the largest. If we 
calculate M e g according to Equation [M] using the best- 
fit parameters in our model, then M c tf = 3x 10 12 at 
z ~ 2, comparable to values found in IViero et alj (l2009f ) 
and IXia et all <|2012f) The o nly direct comparison with 
Mpea^ is IShang et al.1 (|2012f ). Our M pea k is consistent 
with their value log(M pcak /M Q ) = 13.56 ± 0.17 when 
the normalization factor Lq is allowed to vary freely (i.e. 

their Case 4). 

In iCoorav et al.l (|201Q[ ). where a standard halo model 
is developed to fit the angular correlation functions of 
galaxies brighter than 30 mJy, the minimum halo mass 
scale is log(M min /M ) = 12.6, 12.9 and 13.5 at 250, 350 
and 500 u.m, respectively. These values are much higher 
than our best-fit log(M min /M Q ) = 10.9 ± 0.6, which is 
due to the fact that faint galaxies (around 5 mJy) dom- 
inate the power s pectrum of the intens ity fluctuations 
(see Figure [J). In lAmblard et al.l (|201lD . the minimum 
halo mass scale is log(M m j n /MQ) = 11.1, 11.5 and 11.8 at 
250, 350 and 500 u.m, respectively, with which our best- 
fit value is consistent. However, the comparison should 
be made with caution, bearing in mind the differences 



6. DISCUSSION 

We find a clear signature from the clustering of DS- 
FGs in the pattern of CIB anisotropies, with the fidelity 
to identify linear and non-linear terms. Notable is how 
well the spectra are still fit by a power law, with \ 2 ~ 3- 
10 for 15 degrees of freedom (reduced \ 2 ~ 0.3-0.7) — a 
cosmic coincidence first seen in early clustering mea sure- 
ments of resolved galaxies (e.g.. lWatson et aT1l201lt ). Fu- 
ture measurements on larger scales, which should bracket 
the expected peak of the linear, 2-halo spectra, will even- 
tually rule out the power law as a viable fit. That the 
halo model is well motivated regardless is evidenced by 
the change in the power spectrum with changes in the 
level of source masking: while it appears that the 2-halo 
term is negligibly affected, the 1-halo term is significantly 
reduced. We now discuss plausible interpretations of this 
reduction in power. 

6.1. The Reduction of 1-Halo Power with Masking 

Since the 1-halo term originates from multiple galaxies 
occupying the same halo, the reduction of the 1-halo term 
with more aggressive source masking suggests that some 
fraction of the more luminous resolved SPIRE sources 
arc satellites in massive halos or cluster members. This 
is consistent with the interpretation from clustering mea- 
surements of resolved SPIRE sources which claims that 
~ 14% of sources wi th 5 2 50 > 30 mJy appear as satellites 
(jCoorav et al.ll2010D . 

One may wonder if the reduction in power can be solely 
attributed to a reduction in the Poisson level, but this is 
unlikely. We show this by fitting the unmasked 250 u.m 
spectra with a 2-halo, 1-halo, and Poisson term, and 
then fitting the same terms to spectra with sources above 
50 mJy masked. The 2-halo term is fixed in both for both 
levels of masking, while in the latter fit the 1-halo term 
is first fixed and then allowed to float, in order that the 
two fits can be compared. For a floating 1-halo term, 
X 2 = 10.0 for 15 d.o.f. (x 2 oduccd = 0.7), while for a fixed 
1-halo term X 2 = 98.3 for 16 d.o.f. (xLiuced = thus 
ruling out the possibility that the 1-halo term has not 
been affected by masking. 

That some fraction of more luminous sources are found 
in satellite halos is also expect ed from semi-analytic 
models of DSFGs. For example, IConzalez et al.1 (|2011f ) 
predict 38% of all DSFGs (defined there as S > 1 mJy 

at 850 u.m) and 24% of the most luminous sources 
(S 5 mJy at 850 u.m) are satellites. Furthermore, ob- 
servations of individual groups and clusters universally 
find that when star-forming galaxies are present they 
are located on the outskirts of the massive halos (typ- 
ically defined as the vo l ume between R^qq and -R2001 e.g., 
Saintonge et al.l 120081 : iTran et al.1 120091 iBraglia et ail 



201 if) . Even more support for this scenario comes from 



stacking in the submillimeter at positions of brightest 
cluster galaxies (BCGs), which shows a b ump in emis- 
sion at 0.8 Mpc from the central galaxy (|Coppin et al.l 
1 2 1 lh . In light of these observations, the reduction of 
1-halo power with masking is unsurprising. 
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Fig. 13. — Similar to Figurc lTTl but predictions rather than fits are over-plotted on the cross-frequency power spectra with Scut = 300 mJy 
(upper green circles) and S cu t = 50mjy (lower blue circles). Plotted in the top row panels alongside the data are: the Poisson (horizontal 
dashed lines); 2-halo (steep dotted lines); and 1-halo (nearly horizontal dotted lines) terms. The sum of the three terms are shown as solid 
lines. Plotted in the bottom row panels are the predictions minus the data, divided by the data. A perfect fit would have a null value, 
represented here as a horizontal dotted line. 



6.2. Comparison with Planck 

Released to the archive at about the same 
time, initial CIBA powe r sp ectra from Her- 
schel /SPIRE (lAmblard ct all l20l"ll) and Planck /HFI 
(|Planck Collaboration et all 1201 lbl) were found to be 
discrepant by more than 15%. Th e updated, published 
iPlanck Collaboration et all (|2011bfl paper explored this 
difference in detail by comparing to power spectra of 
SPIRE maps that have had no masking applied (rather 
than the published spectra which masked all pixels 
greater than 50mJy). They found that a discrepancy 
still remained, with SPIRE measurements lower by 
factors of - 1.7 and - 1.2 at 857 and 545 GHz (350 and 
550 |J.m) over the angular range 0.02 < fce/arcmin < 0.07 
(400 < i < 1500). If we compare our spectra, having 
only masked extended sources, to the Planck /HFI 
spectra over the same angular range, we find our spectra 
are also low, but by factors of ~ 1.3±0.06 and 1.4±0.06 
over the same range at 350 and 500 ii.m, respectively. 

They next proposed a way to resolve the discrepancy. 
They showed that the spectra of the two groups could 
be brought into agreement if: (i) no Galactic cirrus is re- 
moved from the SPIRE data; and (ii) the beam surface 
area — which they claim to be overestimated by ~ 4 
and 9% at 350 and 500 |J.m, respectively — is corrected. 
Indeed, we also find that the previous cirrus values were 
overestimated; and from a careful estimate of the beam 
area (described in detail in Appendix |B|). we also find 
a correction, though equaling ~ 2 and 8%. These cor- 
rections have brought the overall offset, particularly on 
scales kg & 0.04 arcmin -1 (£ ^ 900) into agreement, how- 
ever, on small scales they have been unable to close the 



gap entirely. The remaining offset on the largest scales 
can be attributed to systematic calibration uncertainties 
of the two instrumen ts: 7% for SPIRE (AppcndixlBl). an d 
7% for Planck /HFI (IPlanck Collaboration ct~aTll2~bllaj ). 
which may be a product of the very different calibration 
strategies of the two instruments; as well as potential 
systematic uncertainties in from cirrus removal. But it 
is on scales kg }t 0.04arcmin _1 {I £ 900) — scales on 

which the contribution from Galactic cirrus is negligible 
— that excess power in the Planck curves are still cither 
in tension (350 |xm) or do not agree (500 p.m). 

It is tempting to try and attribute the remaining small 
angular scale discrepancy to mismatches in the Poisson 
level due to differences in masking, but this solution is 
unlikely. This is because the addition of bright sources 
has a small effect on the Poisson level at 350 urn, and a 
negligible effect at 500 urn (see solid red line in bottom 
panel of Figure^) . And considering that t he comparisons 
were made to SPIRE spectra in which no (jAmblard et alJ 
120111) . or relatively few (~ 1/4 deg 2 ), sources had been 
masked, it is much more likely that this small difference 
in source masking level has a negligible impact on the 
Poisson level. 

A full understanding of the discrepancy for now re- 
mains elusive. That it is mostly the shape of the spec- 
tra on small angular scales that is discrepant, and that 
some models fit to the Planck data were unable to re- 
produ ce the excess of 1-halo power (e.g., IShang et al.l 
l2012t ). points to potential uncertainties in the Planck 
beam measurement. As those were early results, updated 
measurements in the next release from Planck may pro- 
vide the solution. 
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Fig. 14. — Redshift distribution of emission from the best-fit halo 
model. The lAmblard ct al.l 1 120111) points were estimated in a simi- 
lar manner, finding the best-fit solution to their power spect r a with 
a halo model. The model predictions of IBcthcrmin ct al.l (1201 ll ~) 
and Bcthermin (in prep.) are shown as dotted and d ashed lines, re- 
specti vely. Also shown is are model predictions from lValiantc ct al.l 
1120091 ) . The models unanimously anticipate more of a contribution 
from z > 2 than our best-fit finds. 

6.3. Interpreting Halo Model Results 

Our best-fit halo model represents a step forward by 
being the first to fit the power spectra and number 
counts simultaneously. The added complexity introduced 
by this model is motivated by the inability for simpler, 
"standard " halo m odels to simultaneously fit publishe d 
power spectra (e.g.. iPlanck Collaboration et al.l |201 lbf ) ; 
and yet the model remains relatively simple, e.g., it has 
one fewer parameter — two of which wc fix — than 
the standard model presented in lAmblard et al.l (|201lL 
which varied the Poisson level, M m j n , Mi , a, and 4 dS/ dz 
nodes). Despite this, our reduced x 2 of ^ 2 cannot for- 



mally be claimed to be a good fit, and in fact there are 
some obvious problems with this model. 

The first is that the best-estimate of {3 is largely un- 
physi cal (more commonly /3 ~ 2, e.g., iDraine fc Led 
I1984D . Fixing (3 to 1.8 instead dr ives the temperature to 
be 24 K, which is low but not unreasonable, but also leads 
to a s ignificantly worse fit. And although IShang et alj 
(|2012t ) also found their fit of /3 to be in tension, theirs 
was in the opposite direction, with (3 — 2.7. Though 
there is evi dence in the labor atory of a wide range of j3 
values fe.g.. iMcnv ct al.ll2007t ). it is more likely that our 
low /3 value is a symptom of not enough freedom in the 
model SED. 

Also in tension is our effective mass estimate with 
those repor t ed in the literature (|Moster et alj 120101 : 
IShang et all [20121 IWang et all \20lW . Though semi- 
analytic models predict that more than 70% of present 
day bulge-dominated galaxies are descendants of submil- 
limeter galaxies occupying halos of about this mass, our 
host halo mass is much higher than that constrained by 
the correlation function of resolved sources, which is the 
opposite of what is expected. 

Thus, our model is far from the final say in the in- 
terpretation of CIB anisotropies, as it appears that the 
quality of the data demands a model with additional 
levels of sophistication. Future models may potentially 
adopt a set of SEDs for IR emitting galaxies that are ob- 
servationally motivated, and which could include "nor- 
mal" star-forming galaxies, starbursts, and possibly even 
AGN. Furthermore, they could fit not only the three 
SPIRE bands, but also longer wavelength data from e.g., 
ACT, Planck and SPT, as well as at shorter wavelengths 
from e.g., Spitzer and WISE. By combining long wave- 
length multi-band studies of map-based power spectra 
with discrete object correlations at shorter wavelengths, 
we should be able to build a much more complete pic- 
ture of the relationship between stars, star formation, 
and dark matter halos. 

7. SUMMARY 

We have presented the auto- and cross-frequency power 
spectra of cosmic infrared background anisotropies at 
250, 350, and 500 |xm. The background originates from 
all of the dusty star-forming galaxies in the sky; i.e., 
those which are bright and resolved, as well as those too 
faint to be resolved. 

We found an unambiguous signature from the cluster- 
ing of DSFGs in the pattern of the background light and 
showed that it can be decomposed into linear (or 2-halo) 
power from galaxies in separate halos, and non-linear (1- 
halo) power from multiple central and satellite galaxies 
occupying massive halos. We masked resolved sources 
in stages down to 50 mJy and found an expected reduc- 
tion in the level of Poisson noise, as well as a reduction 
in the 1-halo power. We interpreted the reduction in 1- 
halo power as resulting from some fraction of the most 
luminous sources being satellite galaxies. We also mea- 
sured the cross-correlation of the signal between bands 
and found that maps with more aggressive masking to 
be less correlated, as well as hints of a decreasing cor- 
relation with decreasing angular scale; which would be 
indications of decreased correlations between maps for 
higher- z sources. 

We then attempted to interpret the measurement 
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through the framework of a h alo model, building upon 
and extending the formalism of lShang ct al.l (|2012f ) . The 
model was able to fit the power spectra in all three bands 
as well as measured number counts from the literature, 
all simultaneously. We found that, in the framework of 
this model, star formation is most efficient in a range of 
halo masses centered around log(-Mp 0a k/M Q ) ~ 13.9±0.6 
and o\j M ~ 0.3, which is in tension with other estimates 

from the literature. 

Our measurement has limited power to constrain angu- 
lar scales kg £ 0.2arcmin _1 , due partly to the relatively 

small areas of the individual fields, but mostly the re- 
sult of the filtering performed by the SMAP pipeline. 
The situation will improve dramatically with the ar- 
rival of HeLMS which was designed to constrain the 
turnover of the linear term by targeting the largest modes 
in the s ky, as well as fu ture measurements from H- 
ATLAS (|Eales et al.l I2010D and Planck. Add to that 
cross-frequency correlations over the full range of angular 
scales will from e.g., ACT x SPIRE and SPT x SPIRE, 
and even Planck x Planck, which will provide powerful 
new constraints for models of galaxy evolution, the fu- 



ture indeed holds still more breakthroughs 
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APPENDIX 
A. THE UPDATED SMAP PIPELINE 

The reduction and m ap making algorithms used with HcrMES data have evolved since the description presented in 
iLevenson et al.l (|2010f ) . We review the modifications to the SMAP pipeline which lead to the DR1 (first data release) 
HerMES maps in this appendix. These maps are available for download from HcDaM. 5 

Initial processing for the SMAP pipeline uses the Herschel Interactive Processing Environment (HIPE). For HerMES 
DRI, the HCSS / HIPE use r release version 6.0.3, corresponding to continuous integration build 6.0.2055, was used 
(jOtt et al.1 12006L lOttl l2010l ) . including calibration tr ee version spire_cal_6_l. The processing script calls the Spire 
Photometer Interactive Analysis (SPIA: ISchuldl2l)IlT ) version 1.2. 

In summary, the basic pipeline processing steps that are performed by HIPE arc in order: 

1. Signal jump detection. 

2. Common glitch detection. 

3. Sigma-Kappa glitch detection. 

4. Pointing product generation. 

5. Sigma-Kappa glitch repair. 

6. Electronics low-pass filter correction. 

7. Signal linearization and flux calibration. 

8. Bolometer time response correction. 

Details of all these steps, and the implementation of each of the tasks presented below, can be found in the HIPE 
Owner's Guide 6 . Below we detail the tasks called and any changes to their default arguments. 

The initial SMAP processing executes a custom HIPE script that calls the SPIA tasks SPIALevelO_5, SPI- 
aLevelI Repair, spiaLevel2, spiaSaveObs, and spiaSaveMaps2Fits, in that sequence. The arguments to spi- 
aLevelO_5 are the defaults with the following exceptions: 

• waveDeg is set to "Inactive" , switching off the wavelet deglitcher. 

• sigKapDeg is set to "Active" , switching on the sigma-kappa deglitcher. 

• Kappa is set to "4" , meaning that glitches will be detected above 4<r of the timeline noise. 

• LargeGlitchDiscriminatorTimeConstant is set to "4", providing a higher threshold for detecting large glitches. 

Note that the task SPIALevelO_5 provides only detection of jumps and glitches through flags so that the original data 
still can be inspected later in the processing. The flagged glitches are repaired and thermistor timelines with jumps 
are excluded in the task spiaLevelI Repair. 
The parameters in spiaLevelI Repair are default with the following exceptions: 

• extend is set to "Yes" . This will cut off only half of the turn-around datasets after processing one scan, instead 
of the entire turn-around dataset. Because the subsequent scan will keep the other half of the turn-around data, 
the full turn-around dataset remains in the Level 1 data, extending the coverage area. 

• tempDriftCorr is set to "Off" , disabling the temperature correction based on the signals of the thermistor pixels 
on the bolometer arrays. 

The parameters in SPIALevel2 are default with the following exceptions: 



5 http://hedam.oamp.fr/HerMES/ 

6 http : //herschel . esac . esa . int/Docs/DP/HIPE_4 . 2 . 0/hipeowner . pdf 
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• displayMap is set to "No" , preventing the preview images to pop up during processing. 

• makeBrowselmage is set to "No" to prevent generation of browse images irrelevant to this work. 

Note that the HIPE Level 2 maps are not used in the SMAP pipeline. 

The newly-generated Level 1 datasets are then saved in local pools by task SPIASaveObs. These processed 
time streams are then exported to FITS files using the task exportPalToUfDir; the SMAP code reads 
in the time streams at this point. The code base itself is written in the Interactive Data Language (IDL; 
http://www.exelisvis.com/idl/). SMAP first applies a customized set of masks and bad detector lists and ap- 
pends them to the masks carried over from HIPE. These are appended to as required through the following analysis. 

The SPIRE focal planes experience temperature fluctuations which cause the bolometer signals to drift over time. 
These are largely coherent across the focal plane for each array, and can be large (corresponding to as much as 50 Jy 
over 8 h of observation) . The SPIRE focal planes have sensitive thermistor devices that monitor the temperature to 
~ 0.5 u.K at the same sample rate as the detectors; in normal operation, the temperature is stable over 100 s to 2 (J.K, 
so though the instantaneous measurement of the drift is poor, over ~100s scan lengths the signal-to-noise ratio of 
the measurement is > 10. Since the thermistors experience the same fluctuations as the detectors, they can be used 
to remove the component of the bolometer signal arising from the thermal drift in the focal plane. In SMAP, this 
is achieved by stitching together all of the astronomical observation requests (AORs) in a contiguous observation of 
a given field. Both the bolometer and thermistor signals are low-pass filtered with a first-order Butterworth filter 
with a characteristic scale of 1 degree on the sky. Because each SPIRE array has two thermistors, and because the 
thermistors occasionally experience cosmic ray hits or glitches, during times when both thermistors have clean signal 
they are averaged together to improve the signal-to-noise ratio of the measurement of the fluctuations. When one 
thermistor is masked due to data quality issues, the other is used for the duration of the mask. The reconstructed 
average thermistor time stream is then fit to each bolometer in the detector array, and the resulting scaled version of the 
thermistor signal is subtracted. This procedure effectively removes the component of the signal arising from thermal 
fluctuations to the ~ 10~ 4 level. A consequence of this procedure is that the mean (after masking) is subtracted from 
each scan. 

Finally, some scans which pass automated quality masking but which have low-level but visible problems make it 
into the final maps. The maps are inspected, and scans which contribute obviously artificial structure are masked from 
the map making. 

Once the time str eams are completely c onditioned, maps can be constructed. The SMAP map-maker, SHIM, follows 
the presentation in iLevenson et all (|2010h : we summarize here. Our noise model is 

Sdsj = 9dM(x d3j ,y ds j) + pdsj + N dsj , (Al) 

where Sdsj is the signal for detector d, scan s, and time sample j, gd is the detector gain, 7 M(x, y) is the sky brightness 
in pixel (x,y), Ndsj is the instrument noise, and pdsj is an order n polynomial baseline: 



=$>*(**)'■ (A2) 



The parameters a l ds and, optionally, the detector gains, g d , are iteratively fit to the time stream residuals. At each 
iteration i we calculate the residuals: 

R dsj = Sds 3 - [g^P'^Xdsj , y dsj ) + Pdsj] ■ (A3) 

We first fit each of the aj* by minimizing \ 2 — J2j R-dsj with the gd held fixed to g^ 1 , the values calculated from the 
previous iteration. The g d are then fit by minimizing \ 2 — Y2sj ^dsj with the a ds held fixed. On the first iteration, the 
sky is assumed to be 0.0 and the g d are held fixed to 1.0. The sky map M l (x,y) is the weighted mean of all samples 
falling in each pixel: 

E w ds {Sdsj - Pdsj) fad 
M\x,y)= dsie(x ' v) — : , (A4) 



E 



dsje(x,y) 



where the weights w ds are the inverse variance of the timeline residuals, 

-l 

I N 

-jy Z2 ( R dsj) 

3=1 

with TV the number of samples in scan s. 



(A5) 



7 The gd are in fact the deviations from 1.0 of the detector gains already applied by HIPE. 
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The number of iterations and the iteration on which each of the pd s , 9d and w ds are allowed to vary (if any) are all 
specified as inputs to the map maker. For the current data release, DR1, we run for 20 iterations keeping gains fixed 
to 1.0, and allow the weights to deviate from 1.0 starting on the 10 th iteration. 

The SMAP map-maker also performs glitch detection. In addition to the timestream based sigma-kappa glitch 
detection from the HIPE pre-processing mentioned previously, the SMAP map-maker uses an iterative glitch detection 
and removal algorithm based on map information. Taking advantage of the fact that each pixel in the final map is 
sampled by multiple detectors and scans, the SMAP map-maker builds a model of what each detector should see as 
a function of time, including the polynomial baseline. Timestream samples which disagree with this model by more 
than a specified amount (usually 10cr, where a is computed for each timeline after masking) are flagged and removed 
from subsequent map making iterations. This procedure is only activated after a fixed number of iterations (10, by 
default) in order to allow for the values of p to settle, and then is applied for all subsequent iterations. This approach 
is particularly well suited for the HerMES data, which have a large number of scan repeats. 

Finally, we apply an absolute astrometry correction to the maps. This is measured b y stacking preliminary maps 
on Spitzer MIPS 24 urn source s extracted by Vac cari et al. (in prep.) using the SWIRE (jLonsdale et al.ll2003f) MIPS 
24 urn data reduction pipeline (jShupe et alj|2005h . Astrometric registration of MIPS sources was carried out against 
2MASS, returning a mean absolute deviation of the MIPS-2MASS offset of about 0.5"in both RA and Dec in all fields. 
We first make a "quick" map, running the map-maker for only 10 iterations, then make individual maps for each 
AOR using the parameters determined from the quick map. Each AOR map is stacked on the 24 urn catalog and a 
2D Gaussian is fit to the resulting profile. The distance of the center of the fitted Gaussian to the nominal center of 
the image is taken as an absolute shift in the astrometry. These measured offsets are applied to the detector pointing 
solutions in subsequent map-making runs. We note that we have measured the offsets in all three bands independently, 
but find that the measured shifts are consistent between bands, and thus apply the offsets measured at 250 urn, where 
the resolution is highest, to all three bands. The measured shifts are systematic from AOR to AOR, and are generally 
in the range of 1 — 3". 

B. SPIRE MAP CALIBRATION 

Proper calibration of maps is critical for power spectrum measurements, as any systematic offsets are squared in the 
power spectrum. Here we summarize the calibration and color correction procedures; for a more complete description 
see the SPIRE Observers Manual 8 . 

Since Neptune is very bright, relatively compact (angular size ^ 2.5"), and can be seen above instrumental noise in 

the timestreams, SPIRE fluxes are calibrated in the time doma in by fitt i ng the point-spread function (or beam profile) 
to data and setting the peak values to those expected from the iMorenol (|1998| ) model 9 . This measurement is shown to 
be repeatable at the 2% level, and the quoted uncertainty in the Neptune model is 5%, which is conservative and still 
improving. As these are systematic uncertainties, the quoted uncertainty in the calibration is thus 7%. 

In a SPIRE photometer observation, the property of the source that is directly proportional to source power absorbed 
by the bolometer is the integral over the passband of the flux density weighted_by the instrument Relative Spectral 
Response Function (RSRF). Converting from a RSRF- weighted flux density, S§, to a monochromatic flux density 
requires the adoption of a standard frequency for the band and some assumption about the shape of the source 
spectrum. The approach adopted for SPIRE (and PACS) is to assume that the spectrum is a power law across the 
band defined by the flux density at a standard frequency vq, and a spectral index as 

Ss(v) = SsM f - J , (Bl) 

where vq corresponds to frequency equivalent of the nominal SPIRE wavelengths (i.e., 250, 350 and 500 |im), and 
as = — 1, so that the source has a spectrum vS{v) which is flat across the band. The monochromatic flux density at 
frequency vq, which is what is output by HIPE, is then 



y ° J v as °R tYVC {v)dv 



K^typeSs, (B2) 



where "type" refers to point or extended source. For extended sources, the passband is weighted by an additional 
A 7 to account for its width since the beam size increases with increasing wavelength across it. The exact value of 
7 is dependent on the optics of the instrument: though nominally it is expected that the beam area would increase 
as A 2 , in the limit of a very hard taper (or under-illumination) the illumination on the primary is proportional to A, 
and the FWHM on the sky is wavelength independent. The SPIRE taper is slightly wider than a pixel with tophat 
illumination on the primary, meaning that it lies between the two extremes, but closer to nominal. From the optics 
model it is found that 7 = 1.8. The SPIRE photometer pipeline is based on a point source, i.e., -Kpi P = K^-pfaso) = 
[1.0119, 1.0094, 1.0073] at 250, 350, and 500 urn, respectively. 
For extended sources whose true spectra differs from a power law with as = —1, a color correction, Kc.e — 

8 http: //herschel . esac . esa. int/Docs/SPIRE/html/spire_om.html 

9 Tabulations of the Neptune and Uranus brightness temperatures are available from the ESA Herschel Science Centre 
ftp : //ftp . sciops . esa. int/pub/hsc-calibration 
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TABLE 6 

Beam Nominal and Effective Areas. 



Ki^/Kpip, must be applied, where 



Band 


FWHM 


^measured 


Correction 


A eS 


|om 


(arcsec) 


(steradians) 


Factor 


(steradians) 


250 


18.1 


1.039 X 10" 8 


1.013 


1.053 x 10" 8 


350 


25.2 


1.723 x 10~ 8 


1.004 


1.730 x 10~ 8 


500 


36.6 


3.707 x 10~ 8 


0.995 


3.688 x 10" 8 



Ki,E = ^sky(^o) 



J F sky {v)R^{v)dv' 



(B3) 



We would like to color-correct for the case where -FLky^o) is the infrared background, a n extended source which 
FIRAS showed can described by a modified blackbody with T = 18.5 and (3 = 0.65 (jPuget et al.1 Il996h . We 
estimate using the SPIRE passbands additionally weighted by A 2 and a modified black-body SED, finding 

^4,B = [1-0107,1.0022,1.0029] at 250, 350, and 500 Lim, respectively. We check that these corrections are not sen- 
sitive to the approximation made for the FIRAS SED by varying the temperature ±2 K, finding a negligible change 
of ~ ±0.3%. These values compare well with the color corrections for extended sources given in Figure 5.11 of the 
Observers Manual for CIB spectra approximated as power laws across the passbands, with as ~ (0.3,1.1,1.6). In 
summary, the corrections applied to the maps at 250, 350, and 500 u.m are: 

if FIRAS =^C,E = Ki^/Kpip 

= [0.9988,0.9929,0.9957], (B4) 

and hence, negligible. 

Finally, since power spectra are performed on maps in surface brightness units of Jysr -1 , and SPIRE maps are 
natively produced in units of Jy beam -1 , a conversion factor must be applied to the maps, equal to the inverse of the 
solid angle of the beams, 



4eam= / B {0 , <j>)(0 , fldSl, 



(B5) 

where B(6, <f>) is the normalized beam profile, and dfi is the solid angle element in the direction (8, 4>) 

The beam solid angles are measured from SMAP generated maps of Neptune with pixel sizes of 2", normalized by 
the peak value. The area of the SPIRE beams is calculated by summing the Neptune map pixels and multiplying by 
the pixel area, i.e., 4". Next, we address the contamination from background galaxies. We pick a radius, ro, within 
which to integrate, yielding an integral over an area on the map A$. We then pick a second area to be an annular ring 
with r$ < r < hd, where the inner diameter hd = \/To, yielding an area equal to the the inner area. The inner area 
contains the sum of the response to Neptune and the background galaxies, while the outer annulus is just the sum 
of from the galaxies. Assuming the statistics of the background do not change, the outer integral can be subtracted 
from the inner to remove the effect of the background on the beam area. The resulting beam areas given in the 
second column of Table [6] There are systematic uncertainties associated with this calculations. We estimate the beam 
integral by repeating the measurement but varying the values of the input parameters, r and rn> Varying r by ±1' 
from its nominal value results in a fraction of a percent change in the integrals; while varying rjrj by ±10%, which we 
find dominates the error budget, changes the total area by < 1%. Note that these values are specific to SMAP made 
Neptune maps with the same filtering as was used in the maps used in our study, and as such should not be blindly 
adopted for just any SPIRE map. 

Lastly, the beam effective area is corrected for the difference in illumination of the passband due to the relative 
colors of Neptune and the CIB. Both Neptune and the CIB can be described as modified blackbodies, however, the 
temperature of Neptune is ~ 70 K, while the CIB is ~ 18.5 K. To account for this, the beam areas are corrected by the 
ratio of the integrals of the passbands for extended sources weighted by the two SEDs. The resulting correction factors 
are [1.013, 1.004, 0.995] at 250, 350 and 500 urn, respectively. As anticipated, the correction is highest at 250 Lim, where 
the SED of the CIB peaks. We check for potential systematic errors by varying the CIB temperature, and find changes 
to be at the sub-per cent level. The final e ffective beam areas are quoted in the last column of Table |6] Note that the 
beam areas used by lAmblard et al.1 (|2011[) were 1.03, 1.77, and 3.99 x 10 8 steradians. 

C. CONVERSION TO CMB UNITS 

The flux density unit of convention for infrared, (sub)millimeter, and radio astronomers is the Jansky, defined as: 

Jy = 10~ 26 Wm -2 Hz, (CI) 

and is obtained by integrating over the solid angle of the source. For extended sources, the surface brightness is 
described in Jy per unit solid angle, for example, Jysr -1 . Additionally, the power spectrum unit in this convention is 
given in Jy 2 sr _1 . To convert from Jy 2 beam _1 to Jy 2 sr _1 , SPIRE maps must be divided by the area of the beam. 
Beam areas are presented in Table |U For more details see Appendix iBl of this paper, or § 5.2.9 and Table 5.2 of the 
SPIRE Observers Manual 10 . 
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Fig. 15. — C omb ined spectra vs. flux cut of masked sources, plotted as circles with error bars. The best estimates of the cirrus spectra 
in each field (§ 13.31 1 are removed before combining. 

The convention for CMB units is to report a signal as <5Tcmb; *he deviation from the primordial 2.7255 K blackbody. 
To convert from Jysr -1 to <51cmb in /jK, as a function of frequency: 

(G2) 

SB V 2k fkT CMB \ 2 x 2 e x 98.91 Jysr- 1 x 2 e x 

where it=^{-i^) = — jjK — (^Tp (C3) 

and X= kM^ = 56-79 GBz' (C4) 

(|Fixsenll2009f) . Because the SPIRE passbands have widths of - 30% (jGriffin ct al.ll2010D . and because the CMB 
blackbody at these wavelengths is particularly steep (falling exponentially on the Wien side of the 2.7255 K blackbody), 
the integral of 5B U /5T over the bands is weighted towards lower frequencies; an effect that becomes dramatically more 
pronounced at shorter wavelengths. Ultimately, to convert SPIRE maps in Jysr" 1 to u.Kcmb they must be multiplied 
by factors of 3.664 x 1(T 7 , 1.897 x 1CT 8 , and 2.652 x 1(T 10 

To compare the 350 urn band directly to the same band in BLAST requires a slight color correction, as their passbands 
are not quite the same. This correction, from BLAST to SPIRE, is 0.968 in the maps, or 0.937 in the power spectra. 
At 250 and 500 u.m those conversions are respectively 0.994 and 0.996 in the maps, or 0.989 and 0.992 in the power 
spectra, i.e., negligible. 

To compare the 350 urn and 500 u.m bands to the 857 and 545 GHz (or 350 and 550 u.m) Planck bands also requires 
color corrections due to shifts in the band centers. Those conversions, from Planck /HFI to SPIRE at 350 and 500 urn, 
are 0.99 and 1.14 in the maps, or 0.99 and 1.30 in the power spectra, respectively. 

Lastly, the CMB power spectrum is conventionally reported versus multipole £, while in the (sub)millimeter the 
convention is to report it versus angular wavenumber, kg = 1/A, which is also known as er in the literature, and is 
typically expressed in arcmin -1 . In the small-angle approximation the two are related by I = 

D. ALTERNATIVE MASKING SPECTRA 

As described in § 13.11 an alternative set of spectra are calculated on sets of maps which have all been masked 
similarly. I.e., rather than only mask those sources above the flux density cut in that band, all sources identified at 
250 urn are masked at 350 and 500 u.m as well. Note that as before, sources are masked with circles whose sizes are 
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1.1 • FWHM. I.e., the locations of all masked sources are the same band to band, but the sizes of the masks are not. 

The combined spectra are presented in Figure [15] We expect the impact of the alternative masking scheme to be 
most noticeable in the Poisson and 1-halo terms at 350 and 500 Lim, where previously very few sources were masked. 
Poisson level estimates are given in Table [H We find a decline in the Poisson level for the 50mJy flux cut is 5.5 and 
8.8% at 350 and 500 \xm. 

E. TABLES 



TABLE 7 Combined power spectra for all levels of masking. At 

EACH WAVELENGTH, ONLY SOURCES ABOVE THE FLUX CUT ARE MASKED. 



250 x 250 



250 x 350 



Only Extended 
250 X 500 



Sources Masked 
350 x 350 
2.70 ± 1.69) x 10 5 

1.04 ± 0.42) x 10 5 

5.60 ± 1.07) x 10 4 
3.27 ±0.53) x 10 4 
3.34 ±0.37) x 10 4 
2.12 ±0.22) x 10 4 

1.61 ±0.12) x 10 4 
1.23 ±0.06) x 10 4 
1.03 ± 0.04) x 10 4 
8.79 ± 0.32) x 10 3 
8.09 ± 0.27) x 10 3 
7.58 ±0.24) x 10 3 

7.05 ±0.21) x 10 3 
6.72 ±0.20) x 10 3 
6.42 ±0.18) x 10 3 
6.12 ±0.18) x 10 3 
5.84 ±0.24) x 10 3 



350 x 500 
(1.19 ±0.76) x 10 b 
(6.31 ±3.18) x 10 4 
(3.62 ± 0.86) x 10 4 
(1.87 ±0.37) x 10 4 
(1.85 ±0.23) x 10 4 
(1.26 ±0.14) x 10 4 
(9.54 ±0.74) x 10 3 
(7.35 ± 0.38) x 10 3 
(6.15 ±0.26) x 10 3 
(5.26 ±0.18) x 10 3 
(4.88 ±0.15) x 10 3 
(4.53 ±0.13) x 10 3 
(4.20 ±0.12) x 10 3 
(4.03 ±0.11) x 10 3 
(3.81 ±0.12) x 10 3 
(3.62 ±0.18) x 10 3 
(3.39 ± 0.70) x 10 3 



500 x 500 
(7.41 ± 3.87) x 10 4 
(4.45 ± 1.51) x 10 4 
(2.55 ± 0.46) x 10 4 
(1.49 ± 0.22) x 10 4 
(1.23 ± 0.15) x 10 4 
(8.48 ±0.95) x 10 3 
(6.36 ±0.52) x 10 3 
(4.79 ± 0.27) x 10 3 
(3.95 ±0.19) x 10 3 
(3.37 ±0.13) x 10 3 
(3.15 ± 0.11) x 10 3 
(2.90 ±0.09) x 10 3 
(2.68 ± 0.09) x 10 3 
(2.58 ± 0.10) x 10 3 
(2.43 ± 0.13) x 10 3 
(2.33 ± 0.34) x 10 3 
(2.50 ± 1.79) x 10 3 



(5.15 ±3.02) x 10° 
(1.34 ±0.43) x 10 5 
(7.72 ± 1.70) x 10 4 
(5.18 ±0.94) x 10 4 
(4.40 ± 0.54) x 10 4 
(2.84 ±0.32) x 10 4 
(2.27 ±0.20) x 10 4 
(1.67 ±0.11) x 10 4 
(1.45 ±0.08) x 10 4 
(1.23 ±0.06) x 10 4 
(1.13 ±0.05) x 10 4 
(1.08 ± 0.05) x 10 4 
(1.03 ± 0.04) x 10 4 
(9.93 ± 0.39) x 10 3 
(9.40 ± 0.35) x 10 3 
(8.92 ± 0.32) x 10 3 
(8.54 ±0.31) x 10 3 



(3.59 ± 2.45) x 10° 
(1.09 ±0.73) x 10 s 
(5.80 ± 1.57) x 10 4 
(3.39 ± 0.67) x 10 4 
(3.55 ± 0.46) x 10 4 
(2.27 ±0.26) x 10 4 
(1.78 ± 0.14) x 10 4 
(1.35 ± 0.08) x 10 4 
(1.16 ± 0.05) x 10 4 
(9.90 ± 0.42) x 10 3 
(9.10 ± 0.36) x 10 3 
(8.61 ± 0.33) x 10 3 
(8.09 ± 0.28) x 10 3 
(7.77 ±0.26) x 10 3 
(7.38 ±0.23) x 10 3 
(7.01 ± 0.22) x 10 3 
(6.71 ± 0.23) x 10 3 



(1.46 ± 1.08) 
(6.28 ±4.37) 
(4.27 ± 1.08) 
(1.96 ± 0.46) 
(1.95 ±0.28) 
(1.30 ±0.17) 
(9.98 ± 0.86) 
(7.60 ± 0.44) 
(6.47 ±0.31) 
(5.48 ±0.23) 
(5.08 ±0.20) 
(4.74 ±0.17) 
(4.44 ±0.16) 
(4.29 ±0.16) 
(4.03 ±0.17) 
(3.80 ± 0.25) 
(3.59 ±0.86) 



x 10 s 
x 10 4 
x 10 4 
x 10 4 
x 10 4 
x 10 4 
x 10 3 
x 10 3 
x 10 3 
x 10 3 
x 10 3 
x 10 3 
x 10 3 
x 10 3 
x 10 3 
x 10 3 
x 10 3 



250 x 250 



250 x 350 
(3.57 ± 2.44) x 10 b 
(1.09 ± 0.74) x 10 s 
(5.83 ± 1.59) x 10 4 
(3.39 ±0.67) x 10 4 
(3.51 ± 0.46) x 10 4 
(2.20 ±0.26) x 10 4 
(1.73 ± 0.13) x 10 4 
(1.30 ± 0.07) x 10 4 
(1.11 ± 0.04) x 10 4 
(9.48 ± 0.31) x 10 3 
(8.71 ± 0.27) x 10 3 
(8.19 ± 0.24) x 10 3 
(7.69 ± 0.21) x 10 3 
(7.39 ± 0.20) x 10 3 
(7.05 ±0.19) x 10 3 
(6.74 ± 0.18) x 10 3 
(6.50 ± 0.20) x 10 3 



Sources with S > 
250 X 500 



300 mjy Masked 
350 x 350 



350 x 500 



500 x 500 



(5.09 ±2. 
(1.34 ±0. 
(7.73 ± 1. 
(5.11 ±0. 
(4.32 ± 0. 
(2.71 ±0. 
(2.17 ±0. 
(1.57 ±0. 
(1.36 ±0. 
(1.15 ±0. 
(1.06 ±0. 
(1.00 ±0. 
(9.53 ±0. 
(9.22 ±0. 
(8.81 ±0. 
(8.42 ± 0. 
(8.15 ±0. 



99) x 10 J 
43) x 10 5 
70) x 10 4 
93) x 10 4 
52) x 10 4 

30) x 10 4 
17) x 10 4 
08) x 10 4 
06) x 10 4 
04) x 10 4 
04) x 10 4 
04) x 10 4 
32) x 10 3 

31) x 10 3 
29) x 10 3 
27) x 10 3 
27) x 10 3 



(1.45 ± 1.09) 
(6.36 ± 4.41) 
(4.31 ± 1.09) 
(1.95 ±0.45) 
(1.94 ± 0.28) 
(1.26 ±0.16) 
(9.77 ± 0.84) 
(7.42 ± 0.43) 
(6.30 ± 0.28) 
(5.33 ± 0.20) 
(4.92 ± 0.17) 
(4.58 ±0.14) 
(4.29 ±0.13) 
(4.14 ±0.13) 
(3.91 ± 0.14) 
(3.70 ± 0.22) 
(3.51 ± 0.75) 



x 10 s 
x 10 4 
x 10 4 
x 10 4 
x 10 4 
x 10 4 
x 10 3 
x 10 3 
x 10 3 
x 10 3 
x 10 3 
x 10 3 
x 10 3 
x 10 3 
x 10 3 
x 10 3 
x 10 3 



2.70 ± 1.69) 
1.04 ± 0.43) 

5.59 ± 1.07) 
3.29 ±0.53) 
3.32 ±0.38) 
2.11 ±0.23) 

1.60 ±0.12) 
1.22 ±0.06) 
1.02 ± 0.04) 
8.68 ±0.26) 
8.00 ±0.22) 
7.49 ±0.19) 
6.96 ±0.17) 
6.64 ±0.16) 
6.35 ±0.15) 
6.06 ±0.16) 
5.80 ±0.22) 



(1.19 ±0.76) x 10° 

(6.35 ±3.19) x 10 4 

(3.63 ± 0.86) x 10 4 

(1.88 ±0.37) x 10 4 

(1.85 ±0.23) x 10 4 

(1.25 ±0.14) x 10 4 

(9.46 ± 0.73) x 10 3 

(7.31 ±0.39) x 10 3 

(6.10 ±0.25) x 10 3 

(5.21 ±0.17) x 10 3 

(4.84 ±0.14) x 10 3 

x 10 3 (4.49 ±0.12) x 10 3 

(4.16 ± 0.11) x 10 3 

(3.99 ±0.11) x 10 3 

(3.78 ±0.11) x 10 3 

(3.59 ±0.17) x 10 3 

(3.37 ±0.67) x 10 3 



x 10 a 
x 10 5 
X 10 4 
X 10 4 
X 10 4 
X 10 4 
x 10 4 
x 10 4 
x 10 4 
x 10 3 
x 10 3 



x 10 3 

x 

X 
X 
X 



10 3 
10 3 
10 3 
10 3 



(7.41 ± 3.88) x 10 4 
(4.46 ± 1.51) x 10 4 
(2.56 ±0.46) x 10 4 
(1.49 ± 0.22) x 10 4 
(1.23 ± 0.14) x 10 4 
(8.42 ± 0.94) x 10 3 
(6.31 ± 0.51) x 10 3 
(4.77 ±0.28) x 10 3 
(3.93 ±0.18) x 10 3 
(3.35 ± 0.13) x 10 3 
(3.13 ± 0.11) x 10 3 
(2.88 ±0.09) x 10 3 
(2.66 ±0.09) x 10 3 
(2.56 ±0.10) x 10 3 
(2.42 ± 0.13) x 10 3 
(2.32 ±0.34) x 10 3 
(2.49 ± 1.78) x 10 3 



250 x 250 



250 x 350 



Sources with S > 
250 X 500 



200 mjy Masked 
350 X 350 



350 x 500 



500 x 500 



(5.05 ± 2.96) x 10 ;j 
(1.34 ±0.43) x 10 5 
(7.69 ± 1.69) x 10 4 
(5.05 ± 0.92) x 10 4 
(4.26 ± 0.52) x 10 4 
(2.70 ±0.30) x 10 4 
(2.14 ±0.17) x 10 4 
(1.54 ±0.08) x 10 4 
(1.32 ±0.05) x 10 4 
(1.11 ±0.04) x 10 4 
(1.02 ±0.03) x 10 4 
(9.58 ± 0.30) x 10 3 
(9.12 ±0.28) x 10 3 
(8.85 ± 0.27) x 10 3 
(8.46 ± 0.25) x 10 3 
(8.12 ±0.24) x 10 3 
(7.89 ± 0.25) x 10 3 



(3.55 ± 2.43) x 10 s 
(1.09 ± 0.74) x 10 s 
(5.76 ± 1.57) x 10 4 
(3.34 ± 0.66) x 10 4 
(3.48 ± 0.45) x 10 4 
(2.19 ± 0.25) x 10 4 
(1.71 ± 0.13) x 10 4 
(1.29 ±0.06) x 10 4 
(1.10 ± 0.04) x 10 4 
(9.30 ± 0.29) x 10 3 
(8.54 ± 0.24) x 10 3 
(7.99 ± 0.21) x 10 3 
(7.51 ± 0.19) x 10 3 
(7.22 ± 0.18) x 10 3 
(6.90 ±0.17) x 10 3 
(6.60 ± 0.17) x 10 3 
(6.37 ±0.18) x log 



2.70 ± 1.68) x W b 
1.05 ± 0.43) x 10 5 
5.62 ± 1.07) x 10 4 
3.26 ±0.52) x 10 4 
3.31 ±0.37) x 10 4 
2.09 ±0.22) x 10 4 

1.59 ±0.11) x 10 4 
1.21 ±0.06) x 10 4 

1.01 ±0.04) x 10 4 

8.60 ± 0.26) x 10 3 
7.92 ±0.21) x 10 3 
7.41 ±0.18) x 10 3 
6.88 ±0.16) x 10 3 
6.57 ±0.16) x 10 3 
6.30 ±0.15) x 10 3 

6.02 ±0.15) x 10 3 
5.76 ±0.21) x 10 3 



(7.41 ± 3.88) x 10 4 
(4.46 ± 1.51) x 10 4 
(2.56 ±0.46) x 10 4 
(1.49 ±0.22) x 10 4 
(1.23 ± 0.14) x 10 4 
(8.42 ± 0.94) x 10 3 
(6.31 ± 0.52) x 10 3 
(4.77 ±0.28) x 10 3 
(3.93 ± 0.18) x 10 3 
(3.35 ± 0.13) x 10 3 
(3.13 ± 0.11) x 10 3 
(2.88 ± 0.09) x 10 3 
(2.66 ±0.09) x 10 3 
(2.56 ±0.10) x 10 3 
(2.42 ± 0.13) x 10 3 
(2.32 ±0.33) x 10 3 
(2.49 ± 1.77) x 10 3 



(1.44 ± 1.08) x 10 a 
(6.42 ± 4.42) X 10 4 
(4.32 ± 1.09) x 10 4 
(1.91 ± 0.45) x 10 4 
(1.93 ±0.27) x 10 4 
(1.26 ±0.17) x 10 4 
(9.73 ± 0.85) x 10 3 
(7.40 ± 0.42) x 10 3 
(6.23 ± 0.28) x 10 3 
(5.26 ±0.19) x 10 3 
(4.87 ± 0.16) x 10 3 
(4.52 ± 0.13) x 10 3 
(4.23 ±0.12) x 10 3 
(4.08 ±0.12) x 10 3 
(3.85 ±0.13) x 10 3 
(3.65 ± 0.20) x 10 3 
(3.46 ±0.71) x 10 3 



(1.19 ±0.77) x 10" 
(6.36 ± 3.20) x 10 4 
(3.63 ± 0.86) x 10 4 
(1.87 ±0.36) x 10 4 
(1.84 ±0.23) x 10 4 
(1.24 ±0.14) x 10 4 
(9.42 ± 0.74) x 10 3 
(7.27 ±0.39) x 10 3 
(6.08 ±0.25) x 10 3 
(5.19 ±0.17) x 10 3 
(4.81 ±0.14) x 10 3 
(4.46 ±0.12) x 10 3 
(4.13 ±0.10) x 10 3 
(3.96 ±0.10) x 10 3 
(3.75 ±0.11) x 10 3 
(3.57 ±0.16) x 10 3 
(3.36 ±0.65) x 10 3 



Sources with S > 100 mjy Masked 



250 X 250 



250 x 350 



250 x 500 



350 x 350 



350 x 500 



500 x 500 



(4.88 ± 2.91) x W 
(1.33 ± 0.43) x 10 5 
(7.34 ± 1.63) x 10 4 
(4.83 ± 0.89) x 10 4 



(3.44 ± 2.40) x 10 b 
(1.07 ±0.73) x 10 s 
(5.53 ± 1.53) x 10 4 
(3.23 ±0.65) x 10 4 



(2.65 ± 1.67) x W 
(1.04 ±0.42) x 10 5 
(5.46 ± 1.05) x 10 4 
(3.20 ±0.51) x 10 4 



(1.17 ±0.76) x 10 b (7.41 ± 3.88) x 10 4 
(6.37 ± 3.20) x 10 4 (4.45 ± 1.51) x 10 4 
(3.60 ± 0.85) x 10 4 (2.55 ± 0.46) x 10 4 
(1.82 ± 0.36) x 10 4 (1.47 ± 0.22) x 10 4 
Continued on next page. . . 



(1.37 ± 1.06) x 10 J 
(6.31 ±4.39) x 10 4 
(4.31 ± 1.09) x 10 4 
(1.87 ±0.44) x 10 4 
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TABLE 7 - Continued from previous page 



0.044 


(4.20 ± 0.51) 


X 


10 4 


(3.43 ± 0.45) 


X 


10 4 


(1.91 ± 0.27) 


X 


10 4 


(3.28 ± 0.37) 


X 


10 4 


(1.83 ± 


23) 


X 


10 4 


(1.22 ± 0.14) 


X 


10 4 


0.059 


(2.63 ± 0.29) 


X 


10 4 


(2.15 ±0.25) 


X 


10 4 


(1.24 ± 0.16) 


X 


10 4 


(2.06 ± 0.22) 


X 


10 4 


(1.23 ± 


14) 


X 


10 4 


(8.36 ± 0.94) 


X 


10 3 


0.079 


(2.04 ± 0.17) 


X 


10 4 


(1.67 ± 0.13) 


X 


10 4 


(9.54 ± 0.84) 


X 


10 3 


(1.57 ±0.12) 


X 


10 4 


(9.36 ± 


74) 


X 


10 3 


(6.28 ±0.51) 


X 


10 3 


0.105 


(1.46 ± 0.07) 


X 


10 4 


(1.25 ± 0.06) 


X 


10 4 


(7.25 ± 0.41) 


X 


10 3 


(1.19 ± 0.06) 


X 


10 4 


(7.20 ± 


38) 


X 


10 3 


(4.74 ± 0.27) 


X 


10 3 


0.141 


(1.23 ±0.04) 


X 


10 4 


(1.05 ± 0.04) 


X 


10 4 


(6.04 ± 0.27) 


X 


10 3 


(9.84 ± 0.37) 


X 


10 3 


(5.99 ±0 


25) 


X 


10 3 


(3.92 ±0.18) 


X 


10 3 


0.187 


(1.03 ±0.03) 


X 


If) 4 


(8.84 ± 0.26) 


X 


10 3 


(5.10 ±0.18) 


X 


10 3 


(8.35 ± 0.24) 


X 


10 3 


(5.08 ±0 


16) 


X 


10 3 


(3.32 ±0.12) 


X 


10 3 


0.250 


(9.39 ± 0.27) 


X 


10'* 


(8.08 ± 0.21) 


X 


10 3 


(4.69 ±0.14) 


X 


10 3 


(7.68 ±0.19) 


X 


10 3 


(4.71 ± 


13) 


X 


10 3 


(3.11 ± 0.10) 


X 


10 3 


0.333 


(8.79 ± 0.24) 


X 




(7.52 ± 0.18) 


X 


10 3 


(4.34 ±0.12) 


X 


10 3 


(7.14 ±0.17) 


X 


10 3 


(4.35 ± 


11) 


X 


10 3 


(2.85 ±0.09) 


X 


10 3 


n aa^ 


(8.33 ± 0.22) 


X 


i n3 

1U 


(7.05 ± 0.16) 


X 


1U 


(4.05 ±0.11) 


X 


10 3 


(6.65 ±0.15) 


X 


10 3 


(4.03 ± 


10) 


X 


1U 


(2.64 ± 0.08) 


X 


1U 


n ^cn 
u. jyo 


(8.09 ± 0.22) 


X 


i r)3 

±u 


(6.78 ± 0.16) 


X 


1 ri3 


(3.91 ± 0.11) 


X 


10 3 


(6.37 ±0.14) 


X 


10 3 


(3.87 ± 


10) 


X 


i 

1U 


(2.55 ± 0.09) 


X 


1U 


u. / yu 


(7.74 ±0.21) 


X 


i r)3 

1U 


(6.49 ± 0.15) 


X 


1 


(3.69 ±0.12) 


X 


10 3 


(6.11 ±0.14) 


X 


10 3 


(3.68 ± 


10) 


X 


i 

1U 


(2.40 ± 0.13) 


X 


1U 




(7.49 ± 0.20) 


X 


i n3 

1U 


(6.23 ± 0.14) 


X 


1U 


(3.50 ±0.19) 


X 


10 3 


(5.86 ±0.14) 


X 


10 3 


(3.51 ±0 


16) 


X 


1U 


(2.31 ±0.33) 


X 


1U 




(7.31 ±0.21) 


X 


1U 


(6.03 ± 0.16) 


X 


i n3 


(3.31 ± 0.67) 


X 


10 3 


(5.64 ±0.20) 


X 


10 3 


(3.30 ±0 


63) 


X 


i r)3 


(2.48 ± 1.75) 


X 


i 

1U 


kg 














Sources with S > 50mjy Masked 
















[arcmin -1 ] 


250 x 250 




250 x 350 




250 X 500 




350 x 350 




350 x 500 




500 x 500 




0.011 


(4.00 ± 2.62) 


X 


10 5 


(2.55 ± 2.08) 


X 


10 b 


(1.09 ±0.95) 


X 


10 b 


(2.00 ± 1.44) 


X 


10 b 


(9.52 ±6 


76) 


X 


10 4 


(6.80 ± 3.64) 


X 


10 4 


0.019 


(1.25 ±0.41) 


X 


10 5 


(9.88 ± 6.98) 


X 


10 4 


(5.88 ±4.26) 


X 


10 4 


(9.39 ± 3.96) 


X 


10 4 


(5.90 ±3 


03) 


X 


10 4 


(4.22 ± 1.45) 


X 


10 4 


0.026 


(6.64 ± 1.50) 


X 


10 4 


(4.80 ± 1.36) 


X 


10 4 


(4.17 ± 1.05) 


X 


10 4 


(4.94 ± 0.96) 


X 


10 4 


(3.43 ± 


81) 


X 


10 4 


(2.54 ±0.46) 


X 


10 4 


0.033 


(4.46 ± 0.83) 


X 


10 4 


(2.87 ±0.58) 


X 


10 4 


(1.72 ±0.40) 


X 


10 4 


(2.89 ± 0.47) 


X 


10 4 


(1.64 ±0 


33) 


X 


10 4 


(1.38 ±0.21) 


X 


10 4 


0.044 


(3.78 ± 0.46) 


X 


10 4 


(3.08 ± 0.41) 


X 


10 4 


(1.75 ±0.25) 


X 


10 4 


(2.96 ± 0.34) 


X 


10 4 


(1.68 ±0 


21) 


X 


10 4 


(1.17 ±0.14) 


X 


10 4 


0.059 


(2.38 ±0.27) 


X 


10 4 


(1.95 ±0.23) 


X 


10 4 


(1.15 ± 0.15) 


X 


10 4 


(1.88 ±0.20) 


X 


10 4 


(1.13 ±0 


13) 


X 


10 4 


(7.91 ±0.90) 


X 


10 3 


0.079 


(1.77 ±0.14) 


X 


If) 4 


(1.47 ±0.12) 


X 


10 4 


(8.62 ±0.77) 


X 


10 3 


(1.43 ±0.10) 


X 


10 4 


(8.57 ±0 


68) 


X 


10 3 


(5.97 ±0.49) 


X 


10 3 


0.105 


(1.28 ±0.07) 


X 


10 4 


(1.10 ± 0.06) 


X 


10 4 


(6.50 ± 0.38) 


X 


10 3 


(1.07 ±0.05) 


X 


10 4 


(6.50 ± 


34) 


X 


10 3 


(4.43 ± 0.26) 


X 


10 3 


0.141 


(1.03 ± 0.04) 


X 


10 4 


(8.85 ±0.33) 


X 


10 3 


(5.30 ± 0.23) 


X 


10 3 


(8.59 ± 0.32) 


X 


10 3 


(5.31 ±0 


21) 


X 


10 3 


(3.65 ±0.16) 


X 


10 3 


0.187 


(8.57 ±0.26) 


X 


10 3 


(7.48 ±0.21) 


X 


10 3 


(4.44 ±0.15) 


X 


10 3 


(7.40 ±0.21) 


X 


10 3 


(4.53 ± 


14) 


X 


10 3 


(3.11 ±0.11) 


X 


10 3 


0.250 


(7.76 ± 0.22) 


X 


10 3 


(6.79 ± 0.17) 


X 


10 3 


(4.05 ±0.12) 


X 


10 3 


(6.73 ±0.17) 


X 


111 3 


(4.14 ±0 


11) 


X 


10 3 


(2.88 ±0.09) 


X 


10 3 


0.333 


(7.23 ±0.18) 


X 


10 3 


(6.25 ± 0.14) 


X 


10 3 


(3.69 ±0.10) 


X 


10 3 


(6.24 ±0.14) 


X 


111 3 


(3.77 ±0 


09) 


X 


10 3 


(2.63 ± 0.08) 


X 


10 3 


0.445 


(6.76 ±0.17) 


X 


10 3 


(5.79 ± 0.13) 


X 


10 3 


(3.41 ± 0.09) 


X 


10 3 


(5.80 ±0.12) 


X 


111 3 


(3.49 ± 


08) 


X 


10 3 


(2.44 ±0.07) 


X 


10 3 


0.593 


(6.47 ±0.16) 


X 


10 3 


(5.50 ± 0.12) 


X 


10 3 


(3.25 ± 0.09) 


X 


10 3 


(5.50 ±0.12) 


X 


111 3 


(3.32 ± 


08) 


X 


10 3 


(2.36 ± 0.08) 


X 


10 3 


0.790 


(6.18 ±0.15) 


X 


10 3 


(5.25 ± 0.11) 


X 


10 3 


(3.09 ±0.10) 


X 


10 3 


(5.31 ±0.11) 


X 


10 3 


(3.19 ±0 


08) 


X 


10 3 


(2.27 ±0.11) 


X 


10 3 


1.054 


(6.00 ±0.15) 


X 


10 3 


(5.08 ± 0.11) 


X 


10 3 


(2.94 ±0.15) 


X 


10 3 


(5.18 ±0.12) 


X 


10 3 


(3.09 ± 


13) 


X 


10 3 


(2.20 ± 0.29) 


X 


10 3 


1.406 


(5.93 ±0.16) 


X 


10 3 


(5.00 ± 0.13) 


X 


10 3 


(2.80 ±0.53) 


X 


10 3 


(5.05 ±0.17) 


X 


10 3 


(2.96 ±0 


53) 


X 


10 3 


(2.37 ± 1.57) 


X 


10 3 



TABLE 8 

Best-fit Poisson levels as shown in Figure [6] 





S > 50 mjy 


S > 100 mjy 


S > 200 mjy 


S > 300 mjy 


Extended 


250 x 250 
250 x 350 
250 x 500 
350 x 350 
350 x 500 
500 x 500 


(5.97 ±0.05) x 10 3 
(5.06 ± 0.07) x 10 3 
(3.03 ± 0.03) x 10 3 
(5.17 ±0.04) x 10 3 
(3.13 ±0.03) x 10 3 
(2.19 ±0.03) x 10 3 


(7.44 ±0.07) X 10 3 
(6.21 ± 0.09) x 10 3 
(3.64 ±0.04) x 10 3 
(5.90 ±0.09) x 10 3 
(3.63 ± 0.04) x 10 3 
(2.37 ±0.04) x 10 3 


(8.11 ±0.08) x 10 3 
(6.58 ±0.10) x 10 3 
(3.81 ±0.05) x 10 3 
(6.05 ±0.10) x 10 3 
(3.72 ± 0.04) x 10 3 
(2.39 ± 0.04) x 10 3 


(8.42 ± 0.09) X 10 3 
(6.71 ±0.11) x 10 3 
(3.86 ±0.05) x 10 3 
(6.10 ±0.10) x 10 3 
(3.74 ± 0.04) x 10 3 
(2.38 ±0.04) x 10 3 


(8.95 ±0.11) X 10 3 
(6.95 ±0.13) x 10 3 
(3.98 ± 0.06) x 10 3 
(6.16 ±0.12) x 10 3 
(3.77 ± 0.04) x 10 3 
(2.40 ± 0.04) x 10 3 



Note. — At each wavelength, only sources above the flux cut are masked. 



TABLE 9 

Best-fit Poisson levels for spectra with where sources identified at 250 \am are masked in all maps 



S > 50 mjy 



(5.81 ±0.07) x 10 3 " 
(4.98 ±0.05) x 10 3 
(2.90 ±0.06) x 10 3 
(4.67 ±0.05) x 10 3 
(2.90 ±0.05) x 10 3 
(1.94 ±0.05) x 10 3 



S > 100 mjy 



(7.27 ± 0.09) x lO 3 "" 
(6.02 ± 0.07) x 10 3 
(3.45 ± 0.07) x 10 3 
(5.47 ±0.07) x 10 3 
(3.35 ± 0.06) x 10 3 
(2.21 ± 0.06) x 10 3 



S > 200 mjy 



(7.91 ±0.11) X 10 3 
(6.37 ± 0.08) x 10 3 
(3.62 ±0.07) x 10 3 
(5.68 ±0.07) x 10 3 
(3.46 ± 0.06) x 10 3 
(2.26 ±0.06) x 10 3 



S > 300 mjy 



Extended 



(8.20 ±0.12) X lO 3 " 
(6.52 ±0.08) x 10 3 
(3.68 ± 0.08) x 10 3 
(5.75 ±0.07) x 10 3 
(3.49 ± 0.06) x 10 3 
(2.27 ±0.06) x 10 3 



(8.64 ±0.14) x 10 3 
(6.74 ±0.10) x 10 3 
(3.78 ±0.09) x 10 3 
(5.88 ±0.09) x 10 3 
(3.55 ± 0.07) x 10 3 
(2.30 ±0.06) x 10 3 



250 x 250 
250 x 350 
250 x 500 
350 x 350 
350 x 500 
500 x 500 



